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

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