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

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