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

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