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

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