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

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