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

Last change on this file since 5509 was 5505, checked in by gaug, 21 years ago
*** empty log message ***
File size: 32.6 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//
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// Input Containers:
41// MRawEvtData
42// MRawRunHeader
43// MPedestalCam
44//
45// Output Containers:
46// MArrivalTimeCam
47// MExtractedSignalCam
48//
49//////////////////////////////////////////////////////////////////////////////
50#include "MExtractTimeAndChargeDigitalFilter.h"
51
52#include <errno.h>
53#include <fstream>
54
55#include <TFile.h>
56#include <TH1F.h>
57#include <TH2F.h>
58#include <TString.h>
59#include <TMatrix.h>
60
61#include "MLog.h"
62#include "MLogManip.h"
63
64#include "MPedestalPix.h"
65
66ClassImp(MExtractTimeAndChargeDigitalFilter);
67
68using namespace std;
69
70const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst = 0;
71const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast = 14;
72const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst = 3;
73const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast = 14;
74const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeHiGain = 6;
75const Int_t MExtractTimeAndChargeDigitalFilter::fgWindowSizeLoGain = 6;
76const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionHiGain = 10;
77const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolutionLoGain = 10;
78const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinHiGain = 4;
79const Int_t MExtractTimeAndChargeDigitalFilter::fgSignalStartBinLoGain = 4;
80const TString MExtractTimeAndChargeDigitalFilter::fgNameWeightsFile = "cosmics_weights.dat";
81// --------------------------------------------------------------------------
82//
83// Default constructor.
84//
85// Calls:
86// - SetWindowSize();
87// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
88// - SetBinningResolution();
89//
90// Sets all weights to 1.
91//
92MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title)
93 : fWeightsSet(kFALSE), fRandomIter(0)
94{
95 fName = name ? name : "MExtractTimeAndChargeDigitalFilter";
96 fTitle = title ? title : "Digital Filter";
97
98 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
99 SetWindowSize();
100 SetBinningResolution();
101 SetSignalStartBin();
102
103 SetNameWeightsFile();
104}
105
106// ---------------------------------------------------------------------------------------
107//
108// Checks:
109// - if a window is bigger than the one defined by the ranges, set it to the available range
110//
111// Sets:
112// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
113// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
114//
115void MExtractTimeAndChargeDigitalFilter::SetWindowSize(Int_t windowh, Int_t windowl)
116{
117
118 if (windowh != fgWindowSizeHiGain)
119 *fLog << warn << GetDescriptor()
120 << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default." << endl;
121 if (windowl != fgWindowSizeLoGain)
122 *fLog << warn << GetDescriptor()
123 << ": ATTENTION!!! If you are not Hendrik Bartko, do NOT use a different window size than the default" << endl;
124
125 fWindowSizeHiGain = windowh;
126 fWindowSizeLoGain = windowl;
127
128 const Int_t availhirange = (Int_t)(fHiGainLast-fHiGainFirst+1);
129
130 if (fWindowSizeHiGain > availhirange)
131 {
132 // Please simplify this!
133 *fLog << warn << GetDescriptor()
134 << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",fWindowSizeHiGain,
135 " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
136 fHiGainLast = fHiGainFirst + fWindowSizeHiGain;
137 *fLog << warn << GetDescriptor()
138 << ": Will set the upper range to: " << (int)fHiGainLast << endl;
139 }
140
141 if (fWindowSizeHiGain < 2)
142 {
143 fWindowSizeHiGain = 2;
144 *fLog << warn << GetDescriptor() << ": High Gain window size set to two samples" << endl;
145 }
146
147 if (fLoGainLast != 0 && fWindowSizeLoGain != 0)
148 {
149 const Int_t availlorange = (Int_t)(fLoGainLast-fLoGainFirst+1);
150
151 if (fWindowSizeLoGain > availlorange)
152 {
153 // Please simplify this!
154 *fLog << warn << GetDescriptor()
155 << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",fWindowSizeLoGain,
156 " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
157 fLoGainLast = fLoGainFirst + fWindowSizeLoGain;
158 *fLog << warn << GetDescriptor()
159 << ": Will set the upper range to: " << (int)fLoGainLast << endl;
160 }
161
162 if (fWindowSizeLoGain<2)
163 {
164 fWindowSizeLoGain = 2;
165 *fLog << warn << GetDescriptor() << ": Low Gain window size set to two samples" << endl;
166 }
167 }
168
169 //
170 // We need here the effective number of samples which is about 2.5 in the case of a window
171 // size of 6. The exact numbers have to be found still.
172 //
173 fNumHiGainSamples = (Float_t)fWindowSizeHiGain/2.4;
174 fNumLoGainSamples = (Float_t)fWindowSizeLoGain/2.4;
175 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
176 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
177
178}
179
180// --------------------------------------------------------------------------
181//
182// PreProcess
183//
184// Calls:
185// - ReadWeightsFile(fNameWeightsFile) if the weights have not yet been set.
186// - Creates new arrays according to the extraction range
187//
188Int_t MExtractTimeAndChargeDigitalFilter::PreProcess(MParList *pList)
189{
190 if (!fWeightsSet)
191 if (!ReadWeightsFile(fNameWeightsFile))
192 return kFALSE;
193
194 return MExtractTimeAndCharge::PreProcess(pList);
195}
196
197// --------------------------------------------------------------------------
198//
199// ReInit
200//
201// Calls:
202// - MExtractor::ReInit(pList);
203// - Creates new arrays according to the extraction range
204//
205Bool_t MExtractTimeAndChargeDigitalFilter::ReInit(MParList *pList)
206{
207 if (!MExtractTimeAndCharge::ReInit(pList))
208 return kFALSE;
209
210 fTimeShiftHiGain = (Float_t)fHiGainFirst + 0.5 + 1./fBinningResolutionHiGain;
211 fTimeShiftLoGain = (Float_t)fLoGainFirst + 0.5 + 1./fBinningResolutionLoGain;
212
213 return kTRUE;
214}
215
216// --------------------------------------------------------------------------
217//
218// InitArrays
219//
220// Gets called in the ReInit() and initialized the arrays
221//
222void MExtractTimeAndChargeDigitalFilter::InitArrays()
223{
224 Int_t range = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
225 fHiGainSignal.Set(range);
226 range = (Int_t)(fLoGainLast - fLoGainFirst + 1);
227 fLoGainSignal.Set(range);
228}
229
230// --------------------------------------------------------------------------
231//
232// Apply the digital filter algorithm to the high-gain slices.
233//
234void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Float_t &dsum,
235 Float_t &time, Float_t &dtime,
236 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
237{
238
239 Int_t range = fHiGainLast - fHiGainFirst + 1;
240
241 const Byte_t *end = ptr + range;
242 Byte_t *p = ptr;
243 Byte_t maxpos = 0;
244 Byte_t max = 0;
245 Int_t count = 0;
246
247 //
248 // Preparations for the pedestal subtraction (with AB-noise correction)
249 //
250 const Float_t pedes = ped.GetPedestal();
251 const Float_t ABoffs = ped.GetPedestalABoffset();
252
253 Float_t pedmean[2];
254 pedmean[0] = pedes + ABoffs;
255 pedmean[1] = pedes - ABoffs;
256
257 //
258 // Check for saturation in all other slices
259 //
260 while (p<end)
261 {
262
263 const Int_t ids = fHiGainFirst + count ;
264 const Float_t signal = (Float_t)*p - pedmean[(ids+abflag) & 0x1];
265 fHiGainSignal[count] = signal;
266
267 count++;
268
269 if (*p > max)
270 {
271 max = *p;
272 maxpos = p-ptr;
273 }
274
275 if (*p++ >= fSaturationLimit)
276 sat++;
277 }
278
279 if (fHiLoLast != 0)
280 {
281
282 end = logain + fHiLoLast;
283
284 while (logain<end)
285 {
286
287 const Int_t ids = fHiGainFirst + range ;
288 const Float_t signal = (Float_t)*logain - pedmean[(ids+abflag) & 0x1];
289 fHiGainSignal[range] = signal;
290 range++;
291
292 if (*logain > max)
293 {
294 max = *logain;
295 maxpos = range;
296 }
297
298 if (*logain++ >= fSaturationLimit)
299 sat++;
300 }
301 }
302
303 //
304 // allow one saturated slice
305 //
306 if (sat > 0)
307 return;
308
309 //
310 // Slide with a window of size fWindowSizeHiGain over the sample
311 // and multiply the entries with the corresponding weights
312 //
313 if (IsNoiseCalculation())
314 {
315 if (fRandomIter == 10)
316 fRandomIter = 0;
317 for (Int_t ids=0; ids < fWindowSizeHiGain; ids++)
318 {
319 const Int_t idx = fBinningResolutionHiGain*ids + fRandomIter;
320 sum += fAmpWeightsHiGain [idx]*fHiGainSignal[ids];
321 }
322 fRandomIter++;
323 return;
324 }
325
326 Float_t time_sum = 0.;
327 Float_t fmax = 0.;
328 Float_t ftime_max = 0.;
329 Int_t max_p = 0;
330
331 //
332 // Calculate the sum of the first fWindowSize slices
333 //
334 for (Int_t i=0;i<range-fWindowSizeHiGain;i++)
335 {
336 sum = 0.;
337 time_sum = 0.;
338
339 //
340 // Slide with a window of size fWindowSizeHiGain over the sample
341 // and multiply the entries with the corresponding weights
342 //
343 for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
344 {
345 const Int_t idx = fBinningResolutionHiGain*sample+fBinningResolutionHalfHiGain;
346 const Float_t pex = fHiGainSignal[sample+i];
347 sum += fAmpWeightsHiGain [idx]*pex;
348 time_sum += fTimeWeightsHiGain[idx]*pex;
349 }
350
351 if (sum>fmax)
352 {
353 fmax = sum;
354 ftime_max = time_sum;
355 max_p = i;
356 }
357 } /* for (Int_t i=0;i<range-fWindowSizeHiGain;i++) */
358
359 time = 0;
360 if (fmax==0)
361 return;
362
363 ftime_max /= fmax;
364 Int_t t_iter = Int_t(ftime_max*fBinningResolutionHiGain);
365 Int_t sample_iter = 0;
366
367 while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain )
368 {
369 if (t_iter > fBinningResolutionHalfHiGain-1)
370 {
371 t_iter -= fBinningResolutionHiGain;
372 max_p--;
373 sample_iter--;
374 }
375 if (t_iter < -fBinningResolutionHalfHiGain)
376 {
377 t_iter += fBinningResolutionHiGain;
378 max_p++;
379 sample_iter++;
380 }
381 }
382
383 sum = 0.;
384 time_sum = 0.;
385 //
386 // Slide with a window of size fWindowSizeHiGain over the sample
387 // and multiply the entries with the corresponding weights
388 //
389 for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
390 {
391 const Int_t idx = fBinningResolutionHiGain*sample + fBinningResolutionHalfHiGain + t_iter;
392 const Int_t ids = max_p + sample;
393
394 const Float_t pex = ids < 0 ? 0. : ( ids > range ? 0. : fHiGainSignal[ids]);
395 sum += fAmpWeightsHiGain [idx]*pex;
396 time_sum += fTimeWeightsHiGain[idx]*pex;
397 }
398
399 if (sum == 0)
400 return;
401
402 time = max_p + fTimeShiftHiGain /* this shifts the time to the start of the rising edge */
403 - ((Float_t)t_iter)/fBinningResolutionHiGain - time_sum/sum;
404}
405
406// --------------------------------------------------------------------------
407//
408// Apply the digital filter algorithm to the low-gain slices.
409//
410void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
411 Float_t &time, Float_t &dtime,
412 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
413{
414
415 Int_t range = fLoGainLast - fLoGainFirst + 1;
416 const Byte_t *end = ptr + range;
417 Byte_t *p = ptr;
418 Byte_t maxpos = 0;
419 Byte_t max = 0;
420 Int_t count = 0;
421
422 //
423 // Prepare the low-gain pedestal
424 //
425 const Float_t pedes = ped.GetPedestal();
426 const Float_t ABoffs = ped.GetPedestalABoffset();
427
428 Float_t pedmean[2];
429 pedmean[0] = pedes + ABoffs;
430 pedmean[1] = pedes - ABoffs;
431
432 //
433 // Check for saturation in all other slices
434 //
435 while (p<end)
436 {
437
438 const Int_t ids = fLoGainFirst + count ;
439 const Float_t signal = (Float_t)*p - pedmean[(ids+abflag) & 0x1];
440 fLoGainSignal[count] = signal;
441
442 count++;
443
444 if (*p > max)
445 {
446 max = *p;
447 maxpos = p-ptr;
448 }
449
450 if (*p++ >= fSaturationLimit)
451 sat++;
452 }
453
454 //
455 // Slide with a window of size fWindowSizeHiGain over the sample
456 // and multiply the entries with the corresponding weights
457 //
458 if (IsNoiseCalculation())
459 {
460 for (Int_t ids=0; ids < fWindowSizeLoGain; ids++)
461 {
462 const Int_t idx = fBinningResolutionHiGain*ids + fRandomIter;
463 sum += fAmpWeightsLoGain [idx]*fLoGainSignal[ids];
464 }
465 return;
466 }
467
468 Float_t time_sum = 0.;
469 Float_t fmax = 0.;
470 Float_t ftime_max = 0.;
471 Int_t max_p = 0;
472
473 //
474 // Calculate the sum of the first fWindowSize slices
475 //
476 for (Int_t i=0;i<range-fWindowSizeLoGain;i++)
477 {
478 sum = 0.;
479 time_sum = 0.;
480
481 //
482 // Slide with a window of size fWindowSizeLoGain over the sample
483 // and multiply the entries with the corresponding weights
484 //
485 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
486 {
487 const Int_t idx = fBinningResolutionLoGain*sample+fBinningResolutionHalfLoGain;
488 const Float_t pex = fLoGainSignal[sample+i];
489 sum += fAmpWeightsLoGain [idx]*pex;
490 time_sum += fTimeWeightsLoGain[idx]*pex;
491 }
492
493 if (sum>fmax)
494 {
495 fmax = sum;
496 ftime_max = time_sum;
497 max_p = i;
498 }
499 } /* for (Int_t i=0;i<range-fWindowSizeLoGain;i++) */
500
501 time = 0;
502 if (fmax==0)
503 return;
504
505 ftime_max /= fmax;
506 Int_t t_iter = Int_t(ftime_max*fBinningResolutionLoGain);
507 Int_t sample_iter = 0;
508
509 while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain )
510 {
511 if (t_iter > fBinningResolutionHalfLoGain-1)
512 {
513 t_iter -= fBinningResolutionLoGain;
514 max_p--;
515 sample_iter--;
516 }
517 if (t_iter < -fBinningResolutionHalfLoGain)
518 {
519 t_iter += fBinningResolutionLoGain;
520 max_p++;
521 sample_iter++;
522 }
523 }
524
525 sum = 0.;
526 time_sum = 0.;
527
528 //
529 // Slide with a window of size fWindowSizeLoGain over the sample
530 // and multiply the entries with the corresponding weights
531 //
532 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
533 {
534 const Int_t idx = fBinningResolutionLoGain*sample + fBinningResolutionHalfLoGain + t_iter;
535 const Int_t ids = max_p + sample;;
536 const Float_t pex = ids < 0 ? 0. : ( ids > range ? 0. : fLoGainSignal[ids]);
537 sum += fAmpWeightsLoGain [idx]*pex;
538 time_sum += fTimeWeightsLoGain[idx]*pex;
539 }
540
541 if (sum == 0)
542 return;
543
544 time = max_p + fTimeShiftLoGain /* this shifts the time to the start of the rising edge */
545 - ((Float_t)t_iter)/fBinningResolutionLoGain - time_sum/sum;
546}
547
548// --------------------------------------------------------------------------
549//
550// Read the setup from a TEnv, eg:
551// MJPedestal.MExtractor.WindowSizeHiGain: 6
552// MJPedestal.MExtractor.WindowSizeLoGain: 6
553// MJPedestal.MExtractor.BinningResolutionHiGain: 10
554// MJPedestal.MExtractor.BinningResolutionLoGain: 10
555// MJPedestal.MExtractor.WeightsFile: filename
556//
557Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
558{
559 Byte_t hw = fWindowSizeHiGain;
560 Byte_t lw = fWindowSizeLoGain;
561 Bool_t rc = kFALSE;
562
563 if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
564 {
565 hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
566 rc = kTRUE;
567 }
568 if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
569 {
570 lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
571 rc = kTRUE;
572 }
573
574 if (rc)
575 SetWindowSize(hw, lw);
576
577 Bool_t rc2 = kFALSE;
578 Int_t brh = fBinningResolutionHiGain;
579 Int_t brl = fBinningResolutionLoGain;
580
581 if (IsEnvDefined(env, prefix, "BinningResolutionHiGain", print))
582 {
583 brh = GetEnvValue(env, prefix, brh);
584 rc2 = kTRUE;
585 }
586 if (IsEnvDefined(env, prefix, "BinningResolutionLoGain", print))
587 {
588 brl = GetEnvValue(env, prefix, brl);
589 rc2 = kTRUE;
590 }
591
592 if (rc2)
593 {
594 SetBinningResolution(brh, brl);
595 rc = kTRUE;
596 }
597
598 if (IsEnvDefined(env, prefix, "WeightsFile", print))
599 {
600 if (!ReadWeightsFile(GetEnvValue(env, prefix, "WeightsFile", "")))
601 return kERROR;
602 rc = kTRUE;
603 }
604
605 return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
606}
607
608//----------------------------------------------------------------------------
609//
610// Read a pre-defined weights file into the class.
611// This is mandatory for the extraction
612//
613// If filenname is empty, then all weights will be set to 1.
614//
615Bool_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename)
616{
617
618 // This is a fix for TEnv files edited with windows editors
619 filename.ReplaceAll("\015", "");
620
621 SetNameWeightsFile(filename);
622
623 fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain);
624 fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain);
625 fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain);
626 fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain);
627
628 if (fNameWeightsFile.IsNull())
629 {
630 fAmpWeightsHiGain.Reset(1);
631 fTimeWeightsHiGain.Reset(1);
632 fAmpWeightsLoGain.Reset(1);
633 fTimeWeightsLoGain.Reset(1);
634 return kTRUE;
635 }
636
637 ifstream fin(filename.Data());
638 if (!fin)
639 {
640 *fLog << err << GetDescriptor() << ": ERROR - Cannot open file " << filename << ": ";
641 *fLog << strerror(errno) << endl;
642 return kFALSE;
643 }
644
645 *fLog << inf << "Reading weights file " << filename << "..." << flush;
646
647 Int_t len = 0;
648 Int_t cnt = 0;
649 Int_t line = 0;
650 Bool_t hi = kFALSE;
651 Bool_t lo = kFALSE;
652
653 TString str;
654
655 while (1)
656 {
657 str.ReadLine(fin);
658 if (!fin)
659 break;
660
661 line++;
662
663 if (str.Contains("# High Gain Weights:"))
664 {
665 if (hi)
666 {
667 *fLog << err << "ERROR - 'High Gain Weights' found twice in line #" << line << "." << endl;
668 return kFALSE;
669 }
670
671 if (2!=sscanf(str.Data(), "# High Gain Weights:%2i %2i", &fWindowSizeHiGain, &fBinningResolutionHiGain))
672 {
673 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
674 *fLog << str << endl;
675 return kFALSE;
676 }
677
678 len = fBinningResolutionHiGain*fWindowSizeHiGain;
679 fAmpWeightsHiGain .Set(len);
680 fTimeWeightsHiGain.Set(len);
681 hi = kTRUE;
682 continue;
683 }
684
685 if (str.Contains("# Low Gain Weights:"))
686 {
687 if (lo)
688 {
689 *fLog << err << "ERROR - 'Lo Gain Weights' found twice in line #" << line << "." << endl;
690 return kFALSE;
691 }
692
693 if (2!=sscanf(str.Data(),"# Low Gain Weights:%2i %2i", &fWindowSizeLoGain, &fBinningResolutionLoGain))
694 {
695 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
696 *fLog << str << endl;
697 return kFALSE;
698 }
699
700 len = fBinningResolutionLoGain*fWindowSizeHiGain;
701 fAmpWeightsLoGain .Set(len);
702 fTimeWeightsLoGain.Set(len);
703 lo = kTRUE;
704 continue;
705 }
706
707 // Handle lines with comments
708 if (str.Contains("#"))
709 continue;
710
711 // Nothing found so far
712 if (len == 0)
713 continue;
714
715 if (2!=sscanf(str.Data(), "%f %f",
716 lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],
717 lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]))
718 {
719 *fLog << err << "ERROR - Wrong number of arguments in line #" << line << ":" << endl;
720 *fLog << str << endl;
721 return kFALSE;
722 }
723
724 if (++cnt == len)
725 {
726 len = 0;
727 cnt = 0;
728 }
729 }
730
731 if (cnt != len)
732 {
733 *fLog << err << "Size mismatch in weights file " << filename << endl;
734 return kFALSE;
735 }
736
737 if (!hi)
738 {
739 *fLog << err << "No correct header found in weights file " << filename << endl;
740 return kFALSE;
741 }
742
743 *fLog << "done." << endl;
744
745 *fLog << inf << " File contains " << fWindowSizeHiGain << " hi-gain slices ";
746 *fLog << "with a resolution of " << fBinningResolutionHiGain << endl;
747
748 *fLog << inf << " File contains " << fWindowSizeLoGain << " lo-gain slices ";
749 *fLog << "with a resolution of " << fBinningResolutionLoGain << endl;
750
751 fWeightsSet = kTRUE;
752
753 return kTRUE;
754}
755
756//----------------------------------------------------------------------------
757//
758// Create the weights file
759// Beware that the shape-histogram has to contain the pulse starting at bin 1
760//
761Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
762 TH1F *shapelo, TH2F *autocorrlo )
763{
764
765 const Int_t nbinshi = shapehi->GetNbinsX();
766 Float_t binwidth = shapehi->GetBinWidth(1);
767
768 TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
769 Form("%s%s",shapehi->GetTitle()," derivative"),
770 nbinshi,
771 shapehi->GetBinLowEdge(1),
772 shapehi->GetBinLowEdge(nbinshi)+binwidth);
773
774 //
775 // Calculate the derivative of shapehi
776 //
777 for (Int_t i = 1; i<nbinshi+1;i++)
778 {
779 derivativehi->SetBinContent(i,
780 ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
781 derivativehi->SetBinError(i,
782 (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
783 +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
784 }
785
786 //
787 // normalize the shapehi, such that the integral for fWindowSize slices is one!
788 //
789 Float_t sum = 0;
790 Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
791 lasttemp = lasttemp > nbinshi ? nbinshi : lasttemp;
792
793 for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
794 sum += shapehi->GetBinContent(i);
795 }
796 sum /= fBinningResolutionHiGain;
797
798 shapehi->Scale(1./sum);
799 derivativehi->Scale(1./sum);
800
801 //
802 // read in the noise auto-correlation function:
803 //
804 TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);
805
806 for (Int_t i=0; i<fWindowSizeHiGain; i++){
807 for (Int_t j=0; j<fWindowSizeHiGain; j++){
808 Bhi[i][j]=autocorrhi->GetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain
809 }
810 }
811 Bhi.Invert();
812
813 const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
814 fAmpWeightsHiGain.Set(nsizehi);
815 fTimeWeightsHiGain.Set(nsizehi);
816
817 //
818 // Loop over relative time in one BinningResolution interval
819 //
820 Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
821
822 for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++)
823 {
824
825 TMatrix g(fWindowSizeHiGain,1);
826 TMatrix gT(1,fWindowSizeHiGain);
827 TMatrix d(fWindowSizeHiGain,1);
828 TMatrix dT(1,fWindowSizeHiGain);
829
830 for (Int_t count=0; count < fWindowSizeHiGain; count++){
831
832 g[count][0]=shapehi->GetBinContent(start
833 +fBinningResolutionHiGain*count+i);
834 gT[0][count]=shapehi->GetBinContent(start
835 +fBinningResolutionHiGain*count+i);
836 d[count][0]=derivativehi->GetBinContent(start
837 +fBinningResolutionHiGain*count+i);
838 dT[0][count]=derivativehi->GetBinContent(start
839 +fBinningResolutionHiGain*count+i);
840 }
841
842 TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
843 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
844
845 TMatrix m_first = dT*(Bhi*d); // ROOT thinks, m_first is still a matrix
846 Float_t first = m_first[0][0]/denom;
847
848 TMatrix m_last = gT*(Bhi*d); // ROOT thinks, m_last is still a matrix
849 Float_t last = m_last[0][0]/denom;
850
851 TMatrix m1 = gT*Bhi;
852 m1 *= first;
853
854 TMatrix m2 = dT*Bhi;
855 m2 *=last;
856
857 TMatrix w_amp = m1 - m2;
858
859 TMatrix m_first1 = gT*(Bhi*g);
860 Float_t first1 = m_first1[0][0]/denom;
861
862 TMatrix m_last1 = gT*(Bhi*d);
863 Float_t last1 = m_last1 [0][0]/denom;
864
865 TMatrix m11 = dT*Bhi;
866 m11 *=first1;
867
868 TMatrix m21 = gT*Bhi;
869 m21 *=last1;
870
871 TMatrix w_time= m11 - m21;
872
873 for (Int_t count=0; count < fWindowSizeHiGain; count++)
874 {
875 const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1;
876 fAmpWeightsHiGain [idx] = w_amp [0][count];
877 fTimeWeightsHiGain[idx] = w_time[0][count];
878 }
879
880 } // end loop over i
881
882 //
883 // Low Gain histograms
884 //
885 TH1F *derivativelo = NULL;
886 if (shapelo)
887 {
888 const Int_t nbinslo = shapelo->GetNbinsX();
889 binwidth = shapelo->GetBinWidth(1);
890
891 derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
892 Form("%s%s",shapelo->GetTitle()," derivative"),
893 nbinslo,
894 shapelo->GetBinLowEdge(1),
895 shapelo->GetBinLowEdge(nbinslo)+binwidth);
896
897 //
898 // Calculate the derivative of shapelo
899 //
900 for (Int_t i = 1; i<nbinslo+1;i++)
901 {
902 derivativelo->SetBinContent(i,
903 ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
904 derivativelo->SetBinError(i,
905 (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
906 +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
907 }
908
909 //
910 // normalize the shapelo, such that the integral for fWindowSize slices is one!
911 //
912 sum = 0;
913 lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
914 lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
915
916 for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++)
917 sum += shapelo->GetBinContent(i);
918
919 sum /= fBinningResolutionLoGain;
920
921 shapelo->Scale(1./sum);
922 derivativelo->Scale(1./sum);
923
924 //
925 // read in the noise auto-correlation function:
926 //
927 TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
928
929 for (Int_t i=0; i<fWindowSizeLoGain; i++){
930 for (Int_t j=0; j<fWindowSizeLoGain; j++){
931 Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
932 }
933 }
934 Blo.Invert();
935
936 const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
937 fAmpWeightsLoGain.Set(nsizelo);
938 fTimeWeightsLoGain.Set(nsizelo);
939
940 //
941 // Loop over relative time in one BinningResolution interval
942 //
943 Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain;
944
945 for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++)
946 {
947
948 TMatrix g(fWindowSizeLoGain,1);
949 TMatrix gT(1,fWindowSizeLoGain);
950 TMatrix d(fWindowSizeLoGain,1);
951 TMatrix dT(1,fWindowSizeLoGain);
952
953 for (Int_t count=0; count < fWindowSizeLoGain; count++){
954
955 g[count][0] = shapelo->GetBinContent(start
956 +fBinningResolutionLoGain*count+i);
957 gT[0][count]= shapelo->GetBinContent(start
958 +fBinningResolutionLoGain*count+i);
959 d[count][0] = derivativelo->GetBinContent(start
960 +fBinningResolutionLoGain*count+i);
961 dT[0][count]= derivativelo->GetBinContent(start
962 +fBinningResolutionLoGain*count+i);
963 }
964
965 TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
966 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
967
968 TMatrix m_first = dT*(Blo*d); // ROOT thinks, m_first is still a matrix
969 Float_t first = m_first[0][0]/denom;
970
971 TMatrix m_last = gT*(Blo*d); // ROOT thinks, m_last is still a matrix
972 Float_t last = m_last[0][0]/denom;
973
974 TMatrix m1 = gT*Blo;
975 m1 *= first;
976
977 TMatrix m2 = dT*Blo;
978 m2 *=last;
979
980 TMatrix w_amp = m1 - m2;
981
982 TMatrix m_first1 = gT*(Blo*g);
983 Float_t first1 = m_first1[0][0]/denom;
984
985 TMatrix m_last1 = gT*(Blo*d);
986 Float_t last1 = m_last1 [0][0]/denom;
987
988 TMatrix m11 = dT*Blo;
989 m11 *=first1;
990
991 TMatrix m21 = gT*Blo;
992 m21 *=last1;
993
994 TMatrix w_time= m11 - m21;
995
996 for (Int_t count=0; count < fWindowSizeLoGain; count++)
997 {
998 const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1;
999 fAmpWeightsLoGain [idx] = w_amp [0][count];
1000 fTimeWeightsLoGain[idx] = w_time[0][count];
1001 }
1002
1003 } // end loop over i
1004 }
1005
1006 ofstream fn(filename.Data());
1007
1008 fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
1009 fn << "# (Amplitude) (Time) " << endl;
1010
1011 for (Int_t i=0; i<nsizehi; i++)
1012 fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;
1013
1014 fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
1015 fn << "# (Amplitude) (Time) " << endl;
1016
1017 for (Int_t i=0; i<nsizehi; i++)
1018 fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;
1019
1020 delete derivativehi;
1021 if (derivativelo)
1022 delete derivativelo;
1023
1024 return kTRUE;
1025}
1026
1027void MExtractTimeAndChargeDigitalFilter::Print(Option_t *o) const
1028{
1029 *fLog << all;
1030 *fLog << GetDescriptor() << ":" << endl;
1031
1032 MExtractTimeAndCharge::Print(o);
1033 *fLog << " Time Shift HiGain: " << fTimeShiftHiGain << endl;
1034 *fLog << " Time Shift LoGain: " << fTimeShiftLoGain << endl;
1035 *fLog << " Window Size HiGain: " << fWindowSizeHiGain << endl;
1036 *fLog << " Window Size LoGain: " << fWindowSizeLoGain << endl;
1037 *fLog << " Binning Res HiGain: " << fBinningResolutionHiGain << endl;
1038 *fLog << " Binning Res LoGain: " << fBinningResolutionHiGain << endl;
1039
1040 TString opt(o);
1041 if (!opt.Contains("weights"))
1042 return;
1043
1044 *fLog << endl;
1045 *fLog << inf << "Using the following weights: " << endl;
1046 *fLog << "Hi-Gain:" << endl;
1047 for (Int_t i=0; i<fBinningResolutionHiGain*fWindowSizeHiGain; i++)
1048 *fLog << " " << fAmpWeightsHiGain[i] << " \t " << fTimeWeightsHiGain[i] << endl;
1049
1050 *fLog << "Lo-Gain:" << endl;
1051 for (Int_t i=0; i<fBinningResolutionLoGain*fWindowSizeLoGain; i++)
1052 *fLog << " " << fAmpWeightsLoGain[i] << " \t " << fTimeWeightsLoGain[i] << endl;
1053
1054}
Note: See TracBrowser for help on using the repository browser.