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

Last change on this file since 8047 was 7895, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 42.4 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;
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 = -FLT_MAX;
407 Float_t ftime_max = 0.;
408 Int_t max_p = -1;
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 (max_p<0)
439 {
440 sum = 0;
441 time = -1;
442 dsum = -1;
443 dtime = -1;
444 return;
445 }
446
447 if (fmax==0)
448 {
449 sum = 0;
450 time = -1;
451 dsum = 0;
452 dtime = 0;
453 return;
454 }
455
456 ftime_max /= fmax;
457 Int_t t_iter = Int_t(ftime_max*fBinningResolutionHiGain);
458 Int_t sample_iter = 0;
459
460 while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain )
461 {
462 if (t_iter > fBinningResolutionHalfHiGain-1)
463 {
464 t_iter -= fBinningResolutionHiGain;
465 max_p--;
466 sample_iter--;
467 }
468 if (t_iter < -fBinningResolutionHalfHiGain)
469 {
470 t_iter += fBinningResolutionHiGain;
471 max_p++;
472 sample_iter++;
473 }
474 }
475
476 sum = 0.;
477 time_sum = 0.;
478 if (max_p < 0)
479 max_p = 0;
480 if (max_p+fWindowSizeHiGain > range)
481 max_p = range-fWindowSizeHiGain;
482 //
483 // Slide with a window of size fWindowSizeHiGain over the sample
484 // and multiply the entries with the corresponding weights
485 //
486 for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
487 {
488 const Int_t idx = fArrBinningResHalfHiGain[sample] + t_iter;
489 const Int_t ids = max_p + sample;
490 const Float_t pex = fHiGainSignal[ids];
491 sum += fAmpWeightsHiGain [idx]*pex;
492 time_sum += fTimeWeightsHiGain[idx]*pex;
493 }
494
495 if (sum == 0)
496 {
497 time = -1;
498 dsum = 0;
499 dtime = 0;
500 return;
501 }
502
503 // here, the first high-gain slice is already included in the fTimeShiftHiGain
504// time = fTimeShiftHiGain + max_p - Float_t(t_iter)/fBinningResolutionHiGain;
505 time = max_p + fTimeShiftHiGain + (Float_t)fHiGainFirst /* this shifts the time to the start of the rising edge */
506 - ((Float_t)t_iter)/fBinningResolutionHiGain;
507
508 const Float_t timefineadjust = time_sum/sum;
509
510 if (TMath::Abs(timefineadjust) < 4./fBinningResolutionHiGain)
511 time -= timefineadjust;
512
513}
514
515// --------------------------------------------------------------------------
516//
517// Apply the digital filter algorithm to the low-gain slices.
518//
519void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
520 Float_t &time, Float_t &dtime,
521 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
522{
523
524 const Int_t range = fLoGainLast - fLoGainFirst + 1;
525
526 const Byte_t *end = ptr + range;
527 Byte_t *p = ptr;
528 //
529 // Prepare the low-gain pedestal
530 //
531 const Float_t pedes = ped.GetPedestal();
532 const Float_t ABoffs = ped.GetPedestalABoffset();
533
534 const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
535
536 //
537 // Check for saturation in all other slices
538 //
539 Float_t *sample = fLoGainSignal.GetArray();
540 Int_t ids = fLoGainFirst;
541 while (p<end)
542 {
543 *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
544
545 if (*p++ >= fSaturationLimit)
546 sat++;
547 }
548
549 //
550 // Slide with a window of size fWindowSizeLoGain over the sample
551 // and multiply the entries with the corresponding weights
552 //
553 if (IsNoiseCalculation())
554 {
555 if (fRandomIter == fBinningResolutionLoGain)
556 fRandomIter = 0;
557 for (Int_t ids=0; ids < fWindowSizeLoGain; ids++)
558 {
559 const Int_t idx = fArrBinningResLoGain[ids] + fRandomIter;
560 sum += fAmpWeightsLoGain [idx]*fLoGainSignal[ids];
561 }
562 return;
563 }
564
565 Float_t time_sum = 0.;
566 Float_t fmax = -FLT_MAX;
567 Float_t ftime_max = 0.;
568 Int_t max_p = -1;
569
570 //
571 // Calculate the sum of the first fWindowSize slices
572 //
573 for (Int_t i=0;i<range-fWindowSizeLoGain+1;i++)
574 {
575 sum = 0.;
576 time_sum = 0.;
577
578 //
579 // Slide with a window of size fWindowSizeLoGain over the sample
580 // and multiply the entries with the corresponding weights
581 //
582 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
583 {
584 const Int_t idx = fArrBinningResHalfLoGain[sample];
585 const Float_t pex = fLoGainSignal[sample+i];
586 sum += fAmpWeightsLoGain [idx]*pex;
587 time_sum += fTimeWeightsLoGain[idx]*pex;
588 }
589
590 if (sum>fmax)
591 {
592 fmax = sum;
593 ftime_max = time_sum;
594 max_p = i;
595 }
596 } /* for (Int_t i=0;i<range-fWindowSizeLoGain+1;i++) */
597
598 if (max_p<0)
599 {
600 sum = 0;
601 time = -1;
602 dsum = -1;
603 dtime = -1;
604 return;
605 }
606
607 if (fmax==0)
608 {
609 sum = 0;
610 time = -1;
611 dsum = 0;
612 dtime = 0;
613 return;
614 }
615
616 ftime_max /= fmax;
617 Int_t t_iter = Int_t(ftime_max*fBinningResolutionLoGain);
618 Int_t sample_iter = 0;
619
620 while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain )
621 {
622 if (t_iter > fBinningResolutionHalfLoGain-1)
623 {
624 t_iter -= fBinningResolutionLoGain;
625 max_p--;
626 sample_iter--;
627 }
628 if (t_iter < -fBinningResolutionHalfLoGain)
629 {
630 t_iter += fBinningResolutionLoGain;
631 max_p++;
632 sample_iter++;
633 }
634 }
635
636 sum = 0.;
637 time_sum = 0.;
638
639 //
640 // Slide with a window of size fWindowSizeLoGain over the sample
641 // and multiply the entries with the corresponding weights
642 //
643 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
644 {
645 const Int_t idx = fArrBinningResHalfLoGain[sample] + t_iter;
646 const Int_t ids = max_p + sample;
647 const Float_t pex = ids < 0 ? 0. : ( ids >= range ? 0. : fLoGainSignal[ids]);
648 sum += fAmpWeightsLoGain [idx]*pex;
649 time_sum += fTimeWeightsLoGain[idx]*pex;
650 }
651
652 if (sum == 0)
653 {
654 time = -1;
655 dsum = 0;
656 dtime = 0;
657 return;
658 }
659
660 time = max_p + fTimeShiftLoGain + (Float_t)fLoGainFirst /* this shifts the time to the start of the rising edge */
661 - ((Float_t)t_iter)/fBinningResolutionLoGain;
662
663 const Float_t timefineadjust = time_sum/sum;
664
665 if (TMath::Abs(timefineadjust) < 4./fBinningResolutionLoGain)
666 time -= timefineadjust;
667
668}
669
670// --------------------------------------------------------------------------
671//
672// Read the setup from a TEnv, eg:
673// MJPedestal.MExtractor.WindowSizeHiGain: 6
674// MJPedestal.MExtractor.WindowSizeLoGain: 6
675// MJPedestal.MExtractor.BinningResolutionHiGain: 10
676// MJPedestal.MExtractor.BinningResolutionLoGain: 10
677// MJPedestal.MExtractor.WeightsFile: filename
678// MJPedestal.MExtractor.AutomaticWeights: off
679//
680Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
681{
682
683 Byte_t hw = fWindowSizeHiGain;
684 Byte_t lw = fWindowSizeLoGain;
685 Bool_t rc = kFALSE;
686
687 if (IsEnvDefined(env, prefix, "AutomaticWeights", print))
688 {
689 EnableAutomaticWeights(GetEnvValue(env, prefix, "AutomaticWeights", fAutomaticWeights));
690 rc = kTRUE;
691 }
692
693 if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
694 {
695 hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
696 rc = kTRUE;
697 }
698 if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
699 {
700 lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
701 rc = kTRUE;
702 }
703
704 if (rc)
705 SetWindowSize(hw, lw);
706
707 Bool_t rc2 = kFALSE;
708 Int_t brh = fBinningResolutionHiGain;
709 Int_t brl = fBinningResolutionLoGain;
710
711 if (IsEnvDefined(env, prefix, "BinningResolutionHiGain", print))
712 {
713 brh = GetEnvValue(env, prefix, brh);
714 rc2 = kTRUE;
715 }
716 if (IsEnvDefined(env, prefix, "BinningResolutionLoGain", print))
717 {
718 brl = GetEnvValue(env, prefix, brl);
719 rc2 = kTRUE;
720 }
721
722 if (rc2)
723 {
724 SetBinningResolution(brh, brl);
725 rc = kTRUE;
726 }
727
728 if (IsEnvDefined(env, prefix, "WeightsFile", print))
729 {
730 SetNameWeightsFile(GetEnvValue(env, prefix, "WeightsFile", ""));
731 rc = kTRUE;
732 }
733
734 return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
735}
736
737//----------------------------------------------------------------------------
738//
739// If automatic weights are requested, no default weights (name.IsNull())
740// are requested, fRunHeader is available and fRunHeader->IsMonteCarloRun()
741// is true prepend "MC_" in from of the name.
742//
743// return poth+name;
744//
745TString MExtractTimeAndChargeDigitalFilter::CompileWeightFileName(TString path, const TString &name) const
746{
747 if (fAutomaticWeights && !name.IsNull() && fRunHeader && fRunHeader->IsMonteCarloRun())
748 path += "MC_";
749
750 path += name;
751
752 return path;
753}
754
755//----------------------------------------------------------------------------
756//
757// Read a pre-defined weights file into the class.
758// This is mandatory for the extraction
759//
760// If filenname is empty, then all weights will be set to 1.
761//
762// Returns:
763// kTRUE: new weights set
764// kFALSE: no weights set
765// kERROR: error
766//
767Int_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename, TString path)
768{
769 if (filename.IsNull())
770 {
771 fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain);
772 fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain);
773 fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain);
774 fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain);
775
776 fAmpWeightsHiGain.Reset(1);
777 fTimeWeightsHiGain.Reset(1);
778 fAmpWeightsLoGain.Reset(1);
779 fTimeWeightsLoGain.Reset(1);
780 return kTRUE;
781 }
782
783 // Add "MC_" in front of the filename if necessary
784 filename = CompileWeightFileName(path, filename);
785
786 //filename = MJob::ExpandPath(filename);
787
788 if (fNameWeightsFileSet==filename)
789 return kFALSE; // No file read
790
791 ifstream fin(filename.Data());
792 if (!fin)
793 {
794 *fLog << err << GetDescriptor() << ": ERROR - Cannot open file " << filename << ": ";
795 *fLog << strerror(errno) << endl;
796 return kERROR;
797 }
798
799 *fLog << all << GetDescriptor() << ": Reading weights in " << filename << "..." << flush;
800
801 Int_t len = 0;
802 Int_t cnt = 0;
803 Int_t line = 0;
804 Bool_t hi = kFALSE;
805 Bool_t lo = kFALSE;
806
807 TString str;
808
809 while (1)
810 {
811 str.ReadLine(fin);
812 if (!fin)
813 break;
814
815 line++;
816
817 if (str.Contains("# High Gain Weights:"))
818 {
819 if (hi)
820 {
821 *fLog << err << "ERROR - 'High Gain Weights' found twice in line #" << line << "." << endl;
822 return kERROR;
823 }
824
825 if (2!=sscanf(str.Data(), "# High Gain Weights:%2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain))
826 {
827 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
828 *fLog << str << endl;
829 return kERROR;
830 }
831
832 len = fBinningResolutionHiGain*fWindowSizeHiGain;
833 fAmpWeightsHiGain .Set(len);
834 fTimeWeightsHiGain.Set(len);
835 hi = kTRUE;
836 continue;
837 }
838
839 if (str.Contains("# Low Gain Weights:"))
840 {
841 if (lo)
842 {
843 *fLog << err << "ERROR - 'Lo Gain Weights' found twice in line #" << line << "." << endl;
844 return kERROR;
845 }
846
847 if (2!=sscanf(str.Data(),"# Low Gain Weights:%2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain))
848 {
849 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
850 *fLog << str << endl;
851 return kERROR;
852 }
853
854 len = fBinningResolutionLoGain*fWindowSizeLoGain;
855 fAmpWeightsLoGain .Set(len);
856 fTimeWeightsLoGain.Set(len);
857 lo = kTRUE;
858 continue;
859 }
860
861 // Handle lines with comments
862 if (str.Contains("#"))
863 continue;
864
865 // Nothing found so far
866 if (len == 0)
867 continue;
868
869 if (2!=sscanf(str.Data(), "%f %f",
870 lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],
871 lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]))
872 {
873 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
874 *fLog << str << endl;
875 return kERROR;
876 }
877
878 if (++cnt == len)
879 {
880 len = 0;
881 cnt = 0;
882 }
883 }
884
885 if (cnt != len)
886 {
887 *fLog << err << "ERROR - Size mismatch in weights file " << filename << endl;
888 return kERROR;
889 }
890
891 if (!hi)
892 {
893 *fLog << err << "ERROR - No correct header found in weights file " << filename << endl;
894 return kERROR;
895 }
896
897 *fLog << "done." << endl;
898
899 *fLog << inf << " File contains " << fWindowSizeHiGain << " hi-gain slices ";
900 *fLog << "with a resolution of " << fBinningResolutionHiGain << endl;
901
902 *fLog << inf << " File contains " << fWindowSizeLoGain << " lo-gain slices ";
903 *fLog << "with a resolution of " << fBinningResolutionLoGain << endl;
904
905 CalcBinningResArrays();
906
907 switch (fWindowSizeHiGain)
908 {
909 case 4:
910 SetResolutionPerPheHiGain(0.036);
911 break;
912 case 6:
913 SetResolutionPerPheHiGain(0.021);
914 break;
915 default:
916 *fLog << warn << "Could not set the high-gain extractor resolution per phe for window size "
917 << fWindowSizeHiGain << endl;
918 }
919
920 switch (fWindowSizeLoGain)
921 {
922 case 4:
923 SetResolutionPerPheLoGain(0.005);
924 break;
925 case 6:
926 SetResolutionPerPheLoGain(0.004);
927 break;
928 default:
929 *fLog << warn << "Could not set the low-gain extractor resolution per phe for window size "
930 << fWindowSizeLoGain << endl;
931 }
932
933 fNameWeightsFileSet = filename;
934
935 return kTRUE;
936}
937
938
939//----------------------------------------------------------------------------
940//
941// The default (+ prepending possible "MC_") is read for:
942//
943// - RunType: Pedestal (independant of fAutomaticWeights)
944// - fAutomaticWeights disabled
945//
946// if fAutomaticWeights enabled:
947// - fNameWeightsFile.IsNull()
948// - !fCalibPattern
949// - fCalibPattern->GetPulserColor()==MCalibrationCam::kNONE
950//
951// If automatic weights are enabled, the case above didn't take place and
952// fNameWeightsFile starts with "calibration_weights_"
953// - the color (blue, UV) is replaced by the appropriate one
954// taken from the calibration pattern
955//
956// In most cases a debug output is printed. Further output about the color
957// determination can be switched on with debug level > 5;
958//
959// Returns:
960// kFALSE: No new weights set
961// kTRUE: New weights set
962// kERROR: Error
963//
964Int_t MExtractTimeAndChargeDigitalFilter::GetAutomaticWeights()
965{
966 const Ssiz_t pos = fNameWeightsFile.Last('/')+1;
967 const Ssiz_t len = fNameWeightsFile.Length();
968
969 // Split file name in path and name
970 TString path = fNameWeightsFile(0, pos>=0?pos:len);
971 TString name = fNameWeightsFile(pos>=0?pos:0, len);
972
973 // Remove trailing "MC_" for automatic weights
974 if (fAutomaticWeights && name.BeginsWith("MC_"))
975 name.Remove(0, 3);
976
977 // In case of a pedetsal run no calibration pattern can be available
978 // the default weights are always used.
979 if (fRunHeader->GetRunType()==MRawRunHeader::kRTPedestal)
980 {
981 *fLog << dbg << "Pedestal file... using default weights: " << fNameWeightsFile << endl;
982 return ReadWeightsFile(name, path);
983 }
984
985 // If automatic weights are switched off use default weights
986 if (!fAutomaticWeights)
987 {
988 *fLog << dbg << "Automatic weights switched off... using default weights: " << fNameWeightsFile << endl;
989 return ReadWeightsFile(name, path);
990 }
991
992 // If automatic weights are switched on but no filename is given raise error
993 if (fNameWeightsFile.IsNull())
994 {
995 *fLog << err << "ERROR - Cannot get automatic weights without default filename." << endl;
996 return kERROR;
997 }
998
999 // If this is no pedestal run, automatic weights are requested and a
1000 // filename for the weights file is given pedestal-extraction from
1001 // cosmics data is assumed.
1002 if (!fCalibPattern)
1003 {
1004 *fLog << dbg << "No decoded calibration pattern available... using default weights: " << fNameWeightsFile << endl;
1005 return ReadWeightsFile(name, path);
1006 }
1007
1008 const Bool_t debug = gLog.GetDebugLevel()>5;
1009
1010 // If no calibration pattern is available do not change the
1011 // current weighs or current weights file name.
1012 if (fCalibPattern->GetPulserColor()==MCalibrationCam::kNONE)
1013 {
1014 // If we are extracting data and the calibration pattern is kNONE
1015 // we assume that it is a data file without interleaved events
1016 // and calibration pattern information available.
1017 if ((fRunHeader->GetRunType()!=MRawRunHeader::kRTData && !fRunHeader->IsMonteCarloRun()) || debug)
1018 *fLog << dbg << "No calibration color set so far... guessing default: " << fNameWeightsFile << endl;
1019
1020 return ReadWeightsFile(name, path);
1021 }
1022
1023 if (debug)
1024 {
1025 *fLog << dbg << endl;
1026 *fLog << underline << GetDescriptor() << endl;
1027 *fLog << " Trying to get automatic weight for " << fNameWeightsFile << endl;
1028 *fLog << " Run type: ";
1029 }
1030
1031 if (name.BeginsWith("calibration_weights_") && fCalibPattern)
1032 {
1033 if (debug)
1034 *fLog << " Calibration with color " << fCalibPattern->GetPulserColorStr() << ", setting ";
1035 switch (fCalibPattern->GetPulserColor())
1036 {
1037 case MCalibrationCam::kBLUE: // 2
1038 case MCalibrationCam::kGREEN: // 1
1039 if (debug)
1040 *fLog << "blue/green, ";
1041 name.ReplaceAll("UV", "blue");
1042 break;
1043
1044 case MCalibrationCam::kUV: // 3
1045 case MCalibrationCam::kCT1: // 0
1046 if (debug)
1047 *fLog << "UV/CT1, ";
1048 name.ReplaceAll("blue", "UV");
1049 break;
1050 case MCalibrationCam::kNONE:
1051 break;
1052 default: // kNone + etc
1053 *fLog << err << "ERROR - Cannot get automatic weights for " << fCalibPattern->GetPulserColorStr() << endl;
1054 return kERROR;
1055 }
1056 }
1057
1058 return ReadWeightsFile(name, path);
1059}
1060
1061//----------------------------------------------------------------------------
1062//
1063// Create the weights file
1064// Beware that the shape-histogram has to contain the pulse starting at bin 1
1065//
1066Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
1067 TH1F *shapelo, TH2F *autocorrlo )
1068{
1069
1070 const Int_t nbinshi = shapehi->GetNbinsX();
1071 Float_t binwidth = shapehi->GetBinWidth(1);
1072
1073 TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
1074 Form("%s%s",shapehi->GetTitle()," derivative"),
1075 nbinshi,
1076 shapehi->GetBinLowEdge(1),
1077 shapehi->GetBinLowEdge(nbinshi)+binwidth);
1078
1079 //
1080 // Calculate the derivative of shapehi
1081 //
1082 for (Int_t i = 1; i<nbinshi+1;i++)
1083 {
1084 derivativehi->SetBinContent(i,
1085 ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
1086 derivativehi->SetBinError(i,
1087 (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
1088 +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
1089 }
1090
1091 //
1092 // normalize the shapehi, such that the integral for fWindowSize slices is one!
1093 //
1094 Float_t sum = 0;
1095 Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
1096 lasttemp = lasttemp > nbinshi ? nbinshi : lasttemp;
1097
1098 for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
1099 sum += shapehi->GetBinContent(i);
1100 }
1101 sum /= fBinningResolutionHiGain;
1102
1103 shapehi->Scale(1./sum);
1104 derivativehi->Scale(1./sum);
1105
1106 //
1107 // read in the noise auto-correlation function:
1108 //
1109 TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);
1110
1111 for (Int_t i=0; i<fWindowSizeHiGain; i++){
1112 for (Int_t j=0; j<fWindowSizeHiGain; j++){
1113 Bhi[i][j]=autocorrhi->GetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain
1114 }
1115 }
1116 Bhi.Invert();
1117
1118 const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
1119 fAmpWeightsHiGain.Set(nsizehi);
1120 fTimeWeightsHiGain.Set(nsizehi);
1121
1122 //
1123 // Loop over relative time in one BinningResolution interval
1124 //
1125 Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
1126
1127 for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++)
1128 {
1129
1130 TMatrix g(fWindowSizeHiGain,1);
1131 TMatrix gT(1,fWindowSizeHiGain);
1132 TMatrix d(fWindowSizeHiGain,1);
1133 TMatrix dT(1,fWindowSizeHiGain);
1134
1135 for (Int_t count=0; count < fWindowSizeHiGain; count++){
1136
1137 g[count][0]=shapehi->GetBinContent(start
1138 +fBinningResolutionHiGain*count+i);
1139 gT[0][count]=shapehi->GetBinContent(start
1140 +fBinningResolutionHiGain*count+i);
1141 d[count][0]=derivativehi->GetBinContent(start
1142 +fBinningResolutionHiGain*count+i);
1143 dT[0][count]=derivativehi->GetBinContent(start
1144 +fBinningResolutionHiGain*count+i);
1145 }
1146
1147 TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
1148 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
1149
1150 TMatrix m_first = dT*(Bhi*d); // ROOT thinks, m_first is still a matrix
1151 Float_t first = m_first[0][0]/denom;
1152
1153 TMatrix m_last = gT*(Bhi*d); // ROOT thinks, m_last is still a matrix
1154 Float_t last = m_last[0][0]/denom;
1155
1156 TMatrix m1 = gT*Bhi;
1157 m1 *= first;
1158
1159 TMatrix m2 = dT*Bhi;
1160 m2 *=last;
1161
1162 TMatrix w_amp = m1 - m2;
1163
1164 TMatrix m_first1 = gT*(Bhi*g);
1165 Float_t first1 = m_first1[0][0]/denom;
1166
1167 TMatrix m_last1 = gT*(Bhi*d);
1168 Float_t last1 = m_last1 [0][0]/denom;
1169
1170 TMatrix m11 = dT*Bhi;
1171 m11 *=first1;
1172
1173 TMatrix m21 = gT*Bhi;
1174 m21 *=last1;
1175
1176 TMatrix w_time= m11 - m21;
1177
1178 for (Int_t count=0; count < fWindowSizeHiGain; count++)
1179 {
1180 const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1;
1181 fAmpWeightsHiGain [idx] = w_amp [0][count];
1182 fTimeWeightsHiGain[idx] = w_time[0][count];
1183 }
1184
1185 } // end loop over i
1186
1187 //
1188 // Low Gain histograms
1189 //
1190 TH1F *derivativelo = NULL;
1191 if (shapelo)
1192 {
1193 const Int_t nbinslo = shapelo->GetNbinsX();
1194 binwidth = shapelo->GetBinWidth(1);
1195
1196 derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
1197 Form("%s%s",shapelo->GetTitle()," derivative"),
1198 nbinslo,
1199 shapelo->GetBinLowEdge(1),
1200 shapelo->GetBinLowEdge(nbinslo)+binwidth);
1201
1202 //
1203 // Calculate the derivative of shapelo
1204 //
1205 for (Int_t i = 1; i<nbinslo+1;i++)
1206 {
1207 derivativelo->SetBinContent(i,
1208 ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
1209 derivativelo->SetBinError(i,
1210 (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
1211 +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
1212 }
1213
1214 //
1215 // normalize the shapelo, such that the integral for fWindowSize slices is one!
1216 //
1217 sum = 0;
1218 lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
1219 lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
1220
1221 for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++)
1222 sum += shapelo->GetBinContent(i);
1223
1224 sum /= fBinningResolutionLoGain;
1225
1226 shapelo->Scale(1./sum);
1227 derivativelo->Scale(1./sum);
1228
1229 //
1230 // read in the noise auto-correlation function:
1231 //
1232 TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
1233
1234 for (Int_t i=0; i<fWindowSizeLoGain; i++){
1235 for (Int_t j=0; j<fWindowSizeLoGain; j++){
1236 Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
1237 }
1238 }
1239 Blo.Invert();
1240
1241 const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
1242 fAmpWeightsLoGain.Set(nsizelo);
1243 fTimeWeightsLoGain.Set(nsizelo);
1244
1245 //
1246 // Loop over relative time in one BinningResolution interval
1247 //
1248 Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain;
1249
1250 for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++)
1251 {
1252
1253 TMatrix g(fWindowSizeLoGain,1);
1254 TMatrix gT(1,fWindowSizeLoGain);
1255 TMatrix d(fWindowSizeLoGain,1);
1256 TMatrix dT(1,fWindowSizeLoGain);
1257
1258 for (Int_t count=0; count < fWindowSizeLoGain; count++){
1259
1260 g[count][0] = shapelo->GetBinContent(start
1261 +fBinningResolutionLoGain*count+i);
1262 gT[0][count]= shapelo->GetBinContent(start
1263 +fBinningResolutionLoGain*count+i);
1264 d[count][0] = derivativelo->GetBinContent(start
1265 +fBinningResolutionLoGain*count+i);
1266 dT[0][count]= derivativelo->GetBinContent(start
1267 +fBinningResolutionLoGain*count+i);
1268 }
1269
1270 TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
1271 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
1272
1273 TMatrix m_first = dT*(Blo*d); // ROOT thinks, m_first is still a matrix
1274 Float_t first = m_first[0][0]/denom;
1275
1276 TMatrix m_last = gT*(Blo*d); // ROOT thinks, m_last is still a matrix
1277 Float_t last = m_last[0][0]/denom;
1278
1279 TMatrix m1 = gT*Blo;
1280 m1 *= first;
1281
1282 TMatrix m2 = dT*Blo;
1283 m2 *=last;
1284
1285 TMatrix w_amp = m1 - m2;
1286
1287 TMatrix m_first1 = gT*(Blo*g);
1288 Float_t first1 = m_first1[0][0]/denom;
1289
1290 TMatrix m_last1 = gT*(Blo*d);
1291 Float_t last1 = m_last1 [0][0]/denom;
1292
1293 TMatrix m11 = dT*Blo;
1294 m11 *=first1;
1295
1296 TMatrix m21 = gT*Blo;
1297 m21 *=last1;
1298
1299 TMatrix w_time= m11 - m21;
1300
1301 for (Int_t count=0; count < fWindowSizeLoGain; count++)
1302 {
1303 const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1;
1304 fAmpWeightsLoGain [idx] = w_amp [0][count];
1305 fTimeWeightsLoGain[idx] = w_time[0][count];
1306 }
1307
1308 } // end loop over i
1309 }
1310
1311 ofstream fn(filename.Data());
1312
1313 fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
1314 fn << "# (Amplitude) (Time) " << endl;
1315
1316 for (Int_t i=0; i<nsizehi; i++)
1317 fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;
1318
1319 fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
1320 fn << "# (Amplitude) (Time) " << endl;
1321
1322 for (Int_t i=0; i<nsizehi; i++)
1323 fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;
1324
1325 delete derivativehi;
1326 if (derivativelo)
1327 delete derivativelo;
1328
1329 return kTRUE;
1330}
1331
1332void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
1333{
1334 if (IsA()==Class())
1335 *fLog << GetDescriptor() << ":" << endl;
1336
1337 MExtractTimeAndCharge::Print(o);
1338 *fLog << " Time Shift HiGain: " << fTimeShiftHiGain << " LoGain: " << fTimeShiftLoGain << endl;
1339 *fLog << " Window Size HiGain: " << fWindowSizeHiGain << " LoGain: " << fWindowSizeLoGain << endl;
1340 *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << " LoGain: " << fBinningResolutionHiGain << endl;
1341 *fLog << " Weights File desired: " << (fNameWeightsFile.IsNull()?"-":fNameWeightsFile.Data()) << endl;
1342 if (!fNameWeightsFileSet.IsNull())
1343 *fLog << " Weights File set: " << fNameWeightsFileSet << endl;
1344
1345 TString opt(o);
1346 if (!opt.Contains("weights"))
1347 return;
1348
1349 *fLog << endl;
1350 *fLog << inf << "Using the following weights: " << endl;
1351 *fLog << "Hi-Gain:" << endl;
1352 for (Int_t i=0; i<fBinningResolutionHiGain*fWindowSizeHiGain; i++)
1353 *fLog << " " << fAmpWeightsHiGain[i] << " \t " << fTimeWeightsHiGain[i] << endl;
1354
1355 *fLog << "Lo-Gain:" << endl;
1356 for (Int_t i=0; i<fBinningResolutionLoGain*fWindowSizeLoGain; i++)
1357 *fLog << " " << fAmpWeightsLoGain[i] << " \t " << fTimeWeightsLoGain[i] << endl;
1358}
Note: See TracBrowser for help on using the repository browser.