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

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