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

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