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

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