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

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