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

Last change on this file since 5164 was 5164, checked in by hbartko, 20 years ago
*** empty log message ***
File size: 27.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//
29// MExtractTimeAndChargeDigitalFilter
30//
31// Hendrik has promised to write more documentation
32//
33//
34// The following variables have to be set by the derived class and
35// do not have defaults:
36// - fNumHiGainSamples
37// - fNumLoGainSamples
38// - fSqrtHiGainSamples
39// - fSqrtLoGainSamples
40//
41// Input Containers:
42// MRawEvtData
43// MRawRunHeader
44// MPedestalCam
45//
46// Output Containers:
47// MArrivalTimeCam
48// MExtractedSignalCam
49//
50//////////////////////////////////////////////////////////////////////////////
51#include "MExtractTimeAndChargeDigitalFilter.h"
52
53#include "MLog.h"
54#include "MLogManip.h"
55
56#include "MPedestalPix.h"
57
58#include "TFile.h"
59#include "TH1F.h"
60#include "TH2F.h"
61#include "TString.h"
62#include "TMatrix.h"
63
64#include <fstream>
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
94 fName = name ? name : "MExtractTimeAndChargeDigitalFilter";
95 fTitle = title ? title : "Digital Filter";
96
97 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
98 SetWindowSize();
99 SetBinningResolution();
100 SetSignalStartBin();
101
102 ReadWeightsFile("");
103}
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 *fLog << warn << GetDescriptor()
133 << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",fWindowSizeHiGain,
134 " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
135 fHiGainLast = fHiGainFirst + fWindowSizeHiGain;
136 *fLog << warn << GetDescriptor()
137 << ": Will set the upper range to: " << (int)fHiGainLast << endl;
138 }
139
140 if (fWindowSizeHiGain < 2)
141 {
142 fWindowSizeHiGain = 2;
143 *fLog << warn << GetDescriptor() << ": High Gain window size set to two samples" << endl;
144 }
145
146 if (fLoGainLast != 0 && fWindowSizeLoGain != 0)
147 {
148 const Int_t availlorange = (Int_t)(fLoGainLast-fLoGainFirst+1);
149
150 if (fWindowSizeLoGain > availlorange)
151 {
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 fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
168 fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
169
170}
171
172// --------------------------------------------------------------------------
173//
174// ReInit
175//
176// Calls:
177// - MExtractor::ReInit(pList);
178// - Creates new arrays according to the extraction range
179//
180Bool_t MExtractTimeAndChargeDigitalFilter::ReInit(MParList *pList)
181{
182
183 if (!MExtractTimeAndCharge::ReInit(pList))
184 return kFALSE;
185
186 Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
187
188 fHiGainSignal.Set(range);
189
190 range = fLoGainLast - fLoGainFirst + 1;
191
192 fLoGainSignal.Set(range);
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
200
201void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Float_t &dsum,
202 Float_t &time, Float_t &dtime,
203 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
204{
205
206 Int_t range = fHiGainLast - fHiGainFirst + 1;
207 const Byte_t *end = ptr + range;
208 Byte_t *p = ptr;
209 Byte_t maxpos = 0;
210 Byte_t max = 0;
211 Int_t count = 0;
212
213 //
214 // Check for saturation in all other slices
215 //
216 while (p<end)
217 {
218
219 fHiGainSignal[count++] = (Float_t)*p;
220
221 if (*p > max)
222 {
223 max = *p;
224 maxpos = p-ptr;
225 }
226
227 if (*p++ >= fSaturationLimit)
228 sat++;
229 }
230
231 if (fHiLoLast != 0)
232 {
233
234 end = logain + fHiLoLast;
235
236 while (logain<end)
237 {
238
239 fHiGainSignal[count++] = (Float_t)*logain;
240 range++;
241
242 if (*logain > max)
243 {
244 max = *logain;
245 maxpos = range;
246 }
247
248 if (*logain++ >= fSaturationLimit)
249 sat++;
250 }
251 }
252
253 //
254 // allow one saturated slice
255 //
256 if (sat > 0)
257 return;
258
259 Float_t time_sum = 0.;
260 Float_t fmax = 0.;
261 Float_t ftime_max = 0.;
262 Int_t max_p = 0;
263
264 const Float_t pedes = ped.GetPedestal();
265 const Float_t ABoffs = ped.GetPedestalABoffset();
266
267 Float_t PedMean[2];
268 PedMean[0] = pedes + ABoffs;
269 PedMean[1] = pedes - ABoffs;
270
271 //
272 // Calculate the sum of the first fWindowSize slices
273 //
274 for (Int_t i=0;i<range-fWindowSizeHiGain;i++)
275 {
276 sum = 0.;
277 time_sum = 0.;
278
279 //
280 // Slide with a window of size fWindowSizeHiGain over the sample
281 // and multiply the entries with the corresponding weights
282 //
283 for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
284 {
285 const Int_t idx = fBinningResolutionHiGain*sample+fBinningResolutionHalfHiGain;
286 const Float_t pex = fHiGainSignal[sample+i]-PedMean[(sample+i+abflag) & 0x1];
287 sum += fAmpWeightsHiGain [idx]*pex;
288 time_sum += fTimeWeightsHiGain[idx]*pex;
289 }
290
291 if (sum>fmax)
292 {
293 fmax = sum;
294 ftime_max = time_sum;
295 max_p = i;
296 }
297 } /* for (Int_t i=0;i<range-fWindowSizeHiGain;i++) */
298
299 if (fmax!=0)
300 {
301 ftime_max /= fmax;
302 Int_t t_iter = Int_t(ftime_max*fBinningResolutionHiGain);
303 Int_t sample_iter = 0;
304
305 while ( t_iter > fBinningResolutionHalfHiGain-1 || t_iter < -fBinningResolutionHalfHiGain )
306 {
307 if (t_iter > fBinningResolutionHalfHiGain-1)
308 {
309 t_iter -= fBinningResolutionHiGain;
310 max_p--;
311 sample_iter--;
312 }
313 if (t_iter < -fBinningResolutionHalfHiGain)
314 {
315 t_iter += fBinningResolutionHiGain;
316 max_p++;
317 sample_iter++;
318 }
319 }
320
321 sum = 0.;
322 //
323 // Slide with a window of size fWindowSizeHiGain over the sample
324 // and multiply the entries with the corresponding weights
325 //
326 for (Int_t sample=0; sample < fWindowSizeHiGain; sample++)
327 {
328 const Int_t idx = fBinningResolutionHiGain*sample + fBinningResolutionHalfHiGain + t_iter;
329 const Int_t ids = max_p + sample;
330 const Float_t pex = ids < 0 ? 0. :
331 ( ids > range ? 0. : fHiGainSignal[ids]-PedMean[(ids+abflag) & 0x1]);
332 sum += fAmpWeightsHiGain [idx]*pex;
333 time_sum += fTimeWeightsHiGain[idx]*pex;
334 }
335
336 if (sum != 0.)
337 time = max_p + fTimeShiftHiGain /* this shifts the time to the start of the rising edge */
338 - ((Float_t)t_iter)/fBinningResolutionHiGain - time_sum/sum;
339 else
340 time = 0.;
341 } /* if (max!=0) */
342 else
343 time=0.;
344
345 return;
346}
347
348void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum,
349 Float_t &time, Float_t &dtime,
350 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
351{
352
353 Int_t range = fLoGainLast - fLoGainFirst + 1;
354 const Byte_t *end = ptr + range;
355 Byte_t *p = ptr;
356 Byte_t maxpos = 0;
357 Byte_t max = 0;
358 Int_t count = 0;
359
360 //
361 // Check for saturation in all other slices
362 //
363 while (p<end)
364 {
365
366 fLoGainSignal[count++] = (Float_t)*p;
367
368 if (*p > max)
369 {
370 max = *p;
371 maxpos = p-ptr;
372 }
373
374 if (*p++ >= fSaturationLimit)
375 sat++;
376 }
377
378 Float_t time_sum = 0.;
379 Float_t fmax = 0.;
380 Float_t ftime_max = 0.;
381 Int_t max_p = 0;
382
383 const Float_t pedes = ped.GetPedestal();
384 const Float_t ABoffs = ped.GetPedestalABoffset();
385
386 Float_t PedMean[2];
387 PedMean[0] = pedes + ABoffs;
388 PedMean[1] = pedes - ABoffs;
389
390 //
391 // Calculate the sum of the first fWindowSize slices
392 //
393 for (Int_t i=0;i<range-fWindowSizeLoGain;i++)
394 {
395 sum = 0.;
396 time_sum = 0.;
397
398 //
399 // Slide with a window of size fWindowSizeLoGain over the sample
400 // and multiply the entries with the corresponding weights
401 //
402 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
403 {
404 const Int_t idx = fBinningResolutionLoGain*sample+fBinningResolutionHalfLoGain;
405 const Float_t pex = fLoGainSignal[sample+i]-PedMean[(sample+i+abflag) & 0x1];
406 sum += fAmpWeightsLoGain [idx]*pex;
407 time_sum += fTimeWeightsLoGain[idx]*pex;
408 }
409
410 if (sum>fmax)
411 {
412 fmax = sum;
413 ftime_max = time_sum;
414 max_p = i;
415 }
416 } /* for (Int_t i=0;i<range-fWindowSizeLoGain;i++) */
417
418 if (fmax!=0)
419 {
420 ftime_max /= fmax;
421 Int_t t_iter = Int_t(ftime_max*fBinningResolutionLoGain);
422 Int_t sample_iter = 0;
423
424 while ( t_iter > fBinningResolutionHalfLoGain-1 || t_iter < -fBinningResolutionHalfLoGain )
425 {
426 if (t_iter > fBinningResolutionHalfLoGain-1)
427 {
428 t_iter -= fBinningResolutionLoGain;
429 max_p--;
430 sample_iter--;
431 }
432 if (t_iter < -fBinningResolutionHalfLoGain)
433 {
434 t_iter += fBinningResolutionLoGain;
435 max_p++;
436 sample_iter++;
437 }
438 }
439
440 sum = 0.;
441 //
442 // Slide with a window of size fWindowSizeLoGain over the sample
443 // and multiply the entries with the corresponding weights
444 //
445 for (Int_t sample=0; sample < fWindowSizeLoGain; sample++)
446 {
447 const Int_t idx = fBinningResolutionLoGain*sample + fBinningResolutionHalfLoGain + t_iter;
448 const Int_t ids = max_p + sample;
449 const Float_t pex = ids < 0 ? 0. :
450 ( ids > range ? 0. : fLoGainSignal[ids]-PedMean[(ids+abflag) & 0x1]);
451 sum += fAmpWeightsLoGain [idx]*pex;
452 time_sum += fTimeWeightsLoGain[idx]*pex;
453 }
454
455 if (sum != 0.)
456 time = max_p + fTimeShiftLoGain /* this shifts the time to the start of the rising edge */
457 - ((Float_t)t_iter)/fBinningResolutionLoGain - time_sum/sum;
458 else
459 time = 0.;
460 } /* if (max!=0) */
461 else
462 time=0.;
463
464 return;
465}
466
467Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
468{
469
470 Byte_t hw = fWindowSizeHiGain;
471 Byte_t lw = fWindowSizeLoGain;
472 Bool_t rc = kFALSE;
473
474 if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
475 {
476 hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
477 rc = kTRUE;
478 }
479 if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
480 {
481 lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
482 rc = kTRUE;
483 }
484
485 if (rc)
486 SetWindowSize(hw, lw);
487
488 if (IsEnvDefined(env, prefix, "BinningResolution", print))
489 {
490 SetBinningResolution(GetEnvValue(env, prefix, "BinningResolutionHiGain", fBinningResolutionHiGain),
491 GetEnvValue(env, prefix, "BinningResolutionLoGain", fBinningResolutionLoGain));
492 rc = kTRUE;
493 }
494
495 rc = MExtractor::ReadEnv(env, prefix, print) ? kTRUE : rc;
496
497 return rc;
498}
499
500Int_t MExtractTimeAndChargeDigitalFilter::PostProcess()
501{
502
503 *fLog << endl;
504 *fLog << inf << "Used High Gain weights in the extractor: " << endl;
505
506 for (Int_t i=0;i<fBinningResolutionHiGain*fWindowSizeHiGain; i++)
507 *fLog << inf << fAmpWeightsHiGain[i] << " " << fTimeWeightsHiGain[i] << endl;
508
509 *fLog << endl;
510 *fLog << inf << "Used Low Gain weights in the extractor: " << endl;
511 for (Int_t i=0;i<fBinningResolutionLoGain*fWindowSizeLoGain; i++)
512 *fLog << inf << fAmpWeightsLoGain[i] << " " << fTimeWeightsLoGain[i] << endl;
513
514 return kTRUE;
515}
516
517//----------------------------------------------------------------------------
518//
519// Read a pre-defined weights file into the class.
520// This is mandatory for the extraction
521//
522// If filenname is empty, then all weights will be set to 1.
523//
524Bool_t MExtractTimeAndChargeDigitalFilter::ReadWeightsFile(TString filename)
525{
526
527 fAmpWeightsHiGain .Set(fBinningResolutionHiGain*fWindowSizeHiGain);
528 fAmpWeightsLoGain .Set(fBinningResolutionLoGain*fWindowSizeLoGain);
529 fTimeWeightsHiGain.Set(fBinningResolutionHiGain*fWindowSizeHiGain);
530 fTimeWeightsLoGain.Set(fBinningResolutionLoGain*fWindowSizeLoGain);
531
532 if (filename.IsNull())
533 {
534 for (UInt_t i=0; i<fAmpWeightsHiGain.GetSize(); i++)
535 {
536 fAmpWeightsHiGain [i] = 1.;
537 fTimeWeightsHiGain[i] = 1.;
538 }
539 for (UInt_t i=0; i<fAmpWeightsLoGain.GetSize(); i++)
540 {
541 fAmpWeightsLoGain [i] = 1.;
542 fTimeWeightsLoGain[i] = 1.;
543 }
544 return kTRUE;
545 }
546
547 ifstream fin(filename.Data());
548
549 if (!fin)
550 {
551 *fLog << err << GetDescriptor()
552 << ": No weights file found: " << filename << endl;
553 return kFALSE;
554 }
555
556 Int_t len = 0;
557 Int_t cnt = 0;
558 Bool_t hi = kFALSE;
559 Bool_t lo = kFALSE;
560
561 TString str;
562
563 while (1)
564 {
565
566 str.ReadLine(fin);
567 if (!fin)
568 break;
569
570
571 if (str.Contains("# High Gain Weights:"))
572 {
573 str.ReplaceAll("# High Gain Weights:","");
574 sscanf(str.Data(),"%2i%2i",&fWindowSizeHiGain,&fBinningResolutionHiGain);
575 *fLog << inf << "Found number of High Gain slices: " << fWindowSizeHiGain
576 << " and High Gain resolution: " << fBinningResolutionHiGain << endl;
577 len = fBinningResolutionHiGain*fWindowSizeHiGain;
578 fAmpWeightsHiGain .Set(len);
579 fTimeWeightsHiGain.Set(len);
580 hi = kTRUE;
581 continue;
582 }
583
584 if (str.Contains("# Low Gain Weights:"))
585 {
586 str.ReplaceAll("# Low Gain Weights:","");
587 sscanf(str.Data(),"%2i%2i",&fWindowSizeLoGain,&fBinningResolutionLoGain);
588 *fLog << inf << "Found number of Low Gain slices: " << fWindowSizeLoGain
589 << " and Low Gain resolution: " << fBinningResolutionLoGain << endl;
590 len = fBinningResolutionLoGain*fWindowSizeHiGain;
591 fAmpWeightsLoGain .Set(len);
592 fTimeWeightsLoGain.Set(len);
593 lo = kTRUE;
594 continue;
595 }
596
597 if (str.Contains("#"))
598 continue;
599
600 if (len == 0)
601 continue;
602
603 sscanf(str.Data(),"\t%f\t%f",lo ? &fAmpWeightsLoGain [cnt] : &fAmpWeightsHiGain [cnt],
604 lo ? &fTimeWeightsLoGain[cnt] : &fTimeWeightsHiGain[cnt]);
605
606 if (++cnt == len)
607 {
608 len = 0;
609 cnt = 0;
610 }
611 }
612
613 if (cnt != len)
614 {
615 *fLog << err << GetDescriptor()
616 << ": Size mismatch in weights file " << filename << endl;
617 return kFALSE;
618 }
619
620 if (!hi)
621 {
622 *fLog << err << GetDescriptor()
623 << ": No correct header found in weights file " << filename << endl;
624 return kFALSE;
625 }
626
627 return kTRUE;
628}
629
630//----------------------------------------------------------------------------
631//
632// Create the weights file
633// Beware that the shape-histogram has to contain the pulse starting at bin 1
634//
635Bool_t MExtractTimeAndChargeDigitalFilter::WriteWeightsFile(TString filename, TH1F *shapehi, TH2F *autocorrhi,
636 TH1F *shapelo, TH2F *autocorrlo )
637{
638
639 const Int_t nbinshi = shapehi->GetNbinsX();
640 Float_t binwidth = shapehi->GetBinWidth(1);
641
642 TH1F *derivativehi = new TH1F(Form("%s%s",shapehi->GetName(),"_der"),
643 Form("%s%s",shapehi->GetTitle()," derivative"),
644 nbinshi,
645 shapehi->GetBinLowEdge(1),
646 shapehi->GetBinLowEdge(nbinshi)+binwidth);
647
648 //
649 // Calculate the derivative of shapehi
650 //
651 for (Int_t i = 1; i<nbinshi+1;i++)
652 {
653 derivativehi->SetBinContent(i,
654 ((shapehi->GetBinContent(i+1)-shapehi->GetBinContent(i-1))/2./binwidth));
655 derivativehi->SetBinError(i,
656 (sqrt(shapehi->GetBinError(i+1)*shapehi->GetBinError(i+1)
657 +shapehi->GetBinError(i-1)*shapehi->GetBinError(i-1))/2./binwidth));
658 }
659
660 //
661 // normalize the shapehi, such that the integral for fWindowSize slices is one!
662 //
663 Float_t sum = 0;
664 Int_t lasttemp = fBinningResolutionHiGain * (fSignalStartBinHiGain + fWindowSizeHiGain);
665 lasttemp = lasttemp > nbinshi ? nbinshi : lasttemp;
666
667 for (Int_t i=fBinningResolutionHiGain*fSignalStartBinHiGain; i<lasttemp; i++) {
668 sum += shapehi->GetBinContent(i);
669 }
670 sum /= fBinningResolutionHiGain;
671
672 shapehi->Scale(1./sum);
673 derivativehi->Scale(1./sum);
674
675 //
676 // read in the noise auto-correlation function:
677 //
678 TMatrix Bhi(fWindowSizeHiGain,fWindowSizeHiGain);
679
680 for (Int_t i=0; i<fWindowSizeHiGain; i++){
681 for (Int_t j=0; j<fWindowSizeHiGain; j++){
682 Bhi[i][j]=autocorrhi->GetBinContent(i+1,j+1); //+fSignalStartBinHiGain +fSignalStartBinHiGain
683 }
684 }
685 Bhi.Invert();
686
687 const Int_t nsizehi = fWindowSizeHiGain*fBinningResolutionHiGain;
688 fAmpWeightsHiGain.Set(nsizehi);
689 fTimeWeightsHiGain.Set(nsizehi);
690
691 //
692 // Loop over relative time in one BinningResolution interval
693 //
694 Int_t start = fBinningResolutionHiGain*(fSignalStartBinHiGain + 1);
695
696 for (Int_t i = -fBinningResolutionHalfHiGain+1; i<=fBinningResolutionHalfHiGain; i++)
697 {
698
699 TMatrix g(fWindowSizeHiGain,1);
700 TMatrix gT(1,fWindowSizeHiGain);
701 TMatrix d(fWindowSizeHiGain,1);
702 TMatrix dT(1,fWindowSizeHiGain);
703
704 for (Int_t count=0; count < fWindowSizeHiGain; count++){
705
706 g[count][0]=shapehi->GetBinContent(start
707 +fBinningResolutionHiGain*count+i);
708 gT[0][count]=shapehi->GetBinContent(start
709 +fBinningResolutionHiGain*count+i);
710 d[count][0]=derivativehi->GetBinContent(start
711 +fBinningResolutionHiGain*count+i);
712 dT[0][count]=derivativehi->GetBinContent(start
713 +fBinningResolutionHiGain*count+i);
714 }
715
716 TMatrix m_denom = (gT*(Bhi*g))*(dT*(Bhi*d)) - (dT*(Bhi*g))*(dT*(Bhi*g));
717 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
718
719 TMatrix m_first = dT*(Bhi*d); // ROOT thinks, m_first is still a matrix
720 Float_t first = m_first[0][0]/denom;
721
722 TMatrix m_last = gT*(Bhi*d); // ROOT thinks, m_last is still a matrix
723 Float_t last = m_last[0][0]/denom;
724
725 TMatrix m1 = gT*Bhi;
726 m1 *= first;
727
728 TMatrix m2 = dT*Bhi;
729 m2 *=last;
730
731 TMatrix w_amp = m1 - m2;
732
733 TMatrix m_first1 = gT*(Bhi*g);
734 Float_t first1 = m_first1[0][0]/denom;
735
736 TMatrix m_last1 = gT*(Bhi*d);
737 Float_t last1 = m_last1 [0][0]/denom;
738
739 TMatrix m11 = dT*Bhi;
740 m11 *=first1;
741
742 TMatrix m21 = gT*Bhi;
743 m21 *=last1;
744
745 TMatrix w_time= m11 - m21;
746
747 for (Int_t count=0; count < fWindowSizeHiGain; count++)
748 {
749 const Int_t idx = i+fBinningResolutionHalfHiGain+fBinningResolutionHiGain*count-1;
750 fAmpWeightsHiGain [idx] = w_amp [0][count];
751 fTimeWeightsHiGain[idx] = w_time[0][count];
752 }
753
754 } // end loop over i
755
756 //
757 // Low Gain histograms
758 //
759 if (shapelo)
760 {
761 const Int_t nbinslo = shapelo->GetNbinsX();
762 binwidth = shapelo->GetBinWidth(1);
763
764 TH1F *derivativelo = new TH1F(Form("%s%s",shapelo->GetName(),"_der"),
765 Form("%s%s",shapelo->GetTitle()," derivative"),
766 nbinslo,
767 shapelo->GetBinLowEdge(1),
768 shapelo->GetBinLowEdge(nbinslo)+binwidth);
769
770 //
771 // Calculate the derivative of shapelo
772 //
773 for (Int_t i = 1; i<nbinslo+1;i++)
774 {
775 derivativelo->SetBinContent(i,
776 ((shapelo->GetBinContent(i+1)-shapelo->GetBinContent(i-1))/2./binwidth));
777 derivativelo->SetBinError(i,
778 (sqrt(shapelo->GetBinError(i+1)*shapelo->GetBinError(i+1)
779 +shapelo->GetBinError(i-1)*shapelo->GetBinError(i-1))/2./binwidth));
780 }
781
782 //
783 // normalize the shapelo, such that the integral for fWindowSize slices is one!
784 //
785 sum = 0;
786 lasttemp = fBinningResolutionLoGain * (fSignalStartBinLoGain + fWindowSizeLoGain);
787 lasttemp = lasttemp > nbinslo ? nbinslo : lasttemp;
788
789 for (Int_t i=fBinningResolutionLoGain*fSignalStartBinLoGain; i<lasttemp; i++)
790 sum += shapelo->GetBinContent(i);
791
792 sum /= fBinningResolutionLoGain;
793
794 shapelo->Scale(1./sum);
795 derivativelo->Scale(1./sum);
796
797 //
798 // read in the noise auto-correlation function:
799 //
800 TMatrix Blo(fWindowSizeLoGain,fWindowSizeLoGain);
801
802 for (Int_t i=0; i<fWindowSizeLoGain; i++){
803 for (Int_t j=0; j<fWindowSizeLoGain; j++){
804 Blo[i][j]=autocorrlo->GetBinContent(i+1+fSignalStartBinLoGain,j+1+fSignalStartBinLoGain);
805 }
806 }
807 Blo.Invert();
808
809 const Int_t nsizelo = fWindowSizeLoGain*fBinningResolutionLoGain;
810 fAmpWeightsLoGain.Set(nsizelo);
811 fTimeWeightsLoGain.Set(nsizelo);
812
813 //
814 // Loop over relative time in one BinningResolution interval
815 //
816 Int_t start = fBinningResolutionLoGain*fSignalStartBinLoGain + fBinningResolutionHalfLoGain;
817
818 for (Int_t i = -fBinningResolutionHalfLoGain+1; i<=fBinningResolutionHalfLoGain; i++)
819 {
820
821 TMatrix g(fWindowSizeLoGain,1);
822 TMatrix gT(1,fWindowSizeLoGain);
823 TMatrix d(fWindowSizeLoGain,1);
824 TMatrix dT(1,fWindowSizeLoGain);
825
826 for (Int_t count=0; count < fWindowSizeLoGain; count++){
827
828 g[count][0] = shapelo->GetBinContent(start
829 +fBinningResolutionLoGain*count+i);
830 gT[0][count]= shapelo->GetBinContent(start
831 +fBinningResolutionLoGain*count+i);
832 d[count][0] = derivativelo->GetBinContent(start
833 +fBinningResolutionLoGain*count+i);
834 dT[0][count]= derivativelo->GetBinContent(start
835 +fBinningResolutionLoGain*count+i);
836 }
837
838 TMatrix m_denom = (gT*(Blo*g))*(dT*(Blo*d)) - (dT*(Blo*g))*(dT*(Blo*g));
839 Float_t denom = m_denom[0][0]; // ROOT thinks, m_denom is still a matrix
840
841 TMatrix m_first = dT*(Blo*d); // ROOT thinks, m_first is still a matrix
842 Float_t first = m_first[0][0]/denom;
843
844 TMatrix m_last = gT*(Blo*d); // ROOT thinks, m_last is still a matrix
845 Float_t last = m_last[0][0]/denom;
846
847 TMatrix m1 = gT*Blo;
848 m1 *= first;
849
850 TMatrix m2 = dT*Blo;
851 m2 *=last;
852
853 TMatrix w_amp = m1 - m2;
854
855 TMatrix m_first1 = gT*(Blo*g);
856 Float_t first1 = m_first1[0][0]/denom;
857
858 TMatrix m_last1 = gT*(Blo*d);
859 Float_t last1 = m_last1 [0][0]/denom;
860
861 TMatrix m11 = dT*Blo;
862 m11 *=first1;
863
864 TMatrix m21 = gT*Blo;
865 m21 *=last1;
866
867 TMatrix w_time= m11 - m21;
868
869 for (Int_t count=0; count < fWindowSizeLoGain; count++)
870 {
871 const Int_t idx = i+fBinningResolutionHalfLoGain+fBinningResolutionLoGain*count-1;
872 fAmpWeightsLoGain [idx] = w_amp [0][count];
873 fTimeWeightsLoGain[idx] = w_time[0][count];
874 }
875
876 } // end loop over i
877 }
878
879
880 ofstream fn(filename.Data());
881
882 fn << "# High Gain Weights: " << fWindowSizeHiGain << " " << fBinningResolutionHiGain << endl;
883 fn << "# (Amplitude) (Time) " << endl;
884
885 for (Int_t i=0; i<nsizehi; i++)
886 fn << "\t" << fAmpWeightsHiGain[i] << "\t" << fTimeWeightsHiGain[i] << endl;
887
888 fn << "# Low Gain Weights: " << fWindowSizeLoGain << " " << fBinningResolutionLoGain << endl;
889 fn << "# (Amplitude) (Time) " << endl;
890
891 for (Int_t i=0; i<nsizehi; i++)
892 fn << "\t" << fAmpWeightsLoGain[i] << "\t" << fTimeWeightsLoGain[i] << endl;
893
894 return kTRUE;
895}
Note: See TracBrowser for help on using the repository browser.