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

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