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

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