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

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