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

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