1919// ROOT headers
2020#include < TClonesArray.h>
2121#include < TMath.h>
22- #include < array>
2322#include < fstream>
2423#include < iostream>
25- #include < numeric>
26- #include < utility>
2724#include < vector>
2825
2926// FAIR headers
3532#include " R3BActafCalData.h"
3633#include " R3BActafCalPar.h"
3734#include " R3BActafMapped2Cal.h"
35+ #include " R3BActafUtils.h"
3836#include " R3BLogger.h"
3937
40- template <class Cont >
41- int FindMaxPosition (const Cont& signal)
42- {
43- return std::distance (signal.begin (), std::max_element (signal.begin (), signal.end ()));
44- }
45-
46- template <class Cont >
47- double IntegratePulse (const Cont& signal, int maxIdx, double baseline = 0 .)
48- {
49- int left = maxIdx, right = maxIdx;
50- const int size = static_cast <int >(signal.size ());
51-
52- while (left > 0 && signal[left] > baseline)
53- --left;
54-
55- while (right < size - 1 && signal[right] > baseline)
56- ++right;
57-
58- return std::accumulate (signal.begin () + left, signal.begin () + right, 0.0 ) - (right - left) * baseline;
59- }
60-
61- template <class Cont >
62- double ComputeBaselineMean (const Cont& signal, int numBins, bool returnMean = 1 )
63- {
64- numBins = std::min (static_cast <int >(numBins * 0.8 ), static_cast <int >(signal.size ()));
65- if (numBins <= 0 )
66- return 0.0 ;
67-
68- double mean = 0 ;
69-
70- mean = std::accumulate (signal.begin (), signal.begin () + numBins, 0.0 ) / numBins;
71-
72- if (returnMean)
73- return mean;
74-
75- else
76- {
77- double variance = 0.0 ;
78- for (int i = 0 ; i < numBins; ++i)
79- {
80- double diff = signal[i] - mean;
81- variance += diff * diff;
82- }
83- variance /= numBins;
84-
85- double stddev = std::sqrt (variance);
86- return stddev;
87- }
88- }
89-
9038// R3BActafMapped2Cal::Default Constructor --------------------------
9139R3BActafMapped2Cal::R3BActafMapped2Cal ()
9240 : R3BActafMapped2Cal(" R3BActafMapped2Cal" , 1 )
@@ -158,19 +106,11 @@ InitStatus R3BActafMapped2Cal::Init()
158106{
159107 R3BLOG (info, " " );
160108 auto * mgr = FairRootManager::Instance ();
161- if (!mgr)
162- {
163- R3BLOG (fatal, " FairRootManager not found" );
164- return kFATAL ;
165- }
109+ R3BLOG_IF (fatal, !mgr, " FairRootManager not found" );
166110
167111 // INPUT DATA
168112 fActafMappedData = dynamic_cast <TClonesArray*>(mgr->GetObject (" ActafMappedData" ));
169- if (!fActafMappedData )
170- {
171- R3BLOG (fatal, " ActafMappedData not found" );
172- return kFATAL ;
173- }
113+ R3BLOG_IF (fatal, !fActafMappedData , " ActafMappedData not found" );
174114
175115 // OUTPUT DATA
176116 fActafCalData = new TClonesArray (" R3BActafCalData" );
@@ -190,32 +130,6 @@ InitStatus R3BActafMapped2Cal::ReInit()
190130 return kSUCCESS ;
191131}
192132
193- void R3BActafMapped2Cal::ApplySGFilter (std::array<double , ACTAF_BINS >& signal, std::vector<double > coeffs)
194- {
195- auto n = signal.size (), m = coeffs.size ();
196- int half = m / 2 ;
197-
198- std::vector<double > output (n), ext (n + 2 * half);
199-
200- for (int i = 0 ; i < half; i++)
201- ext[i] = signal[0 ];
202- for (int i = 0 ; i < n; i++)
203- ext[i + half] = signal[i];
204- for (int i = 0 ; i < half; i++)
205- ext[n + half + i] = signal[n - 1 ];
206-
207- for (int i = 0 ; i < n; i++)
208- {
209- double sum = 0.0 ;
210- for (int j = 0 ; j < m; j++)
211- sum += coeffs[j] * ext[i + j];
212- output[i] = sum;
213- }
214-
215- for (int i = 0 ; i < n; i++)
216- signal[i] = output[i];
217- }
218-
219133// ----- Public method Execution --------------------------------------------
220134void R3BActafMapped2Cal::Exec (Option_t*)
221135{
@@ -246,44 +160,45 @@ void R3BActafMapped2Cal::Exec(Option_t*)
246160 continue ;
247161
248162 std::array<double , ACTAF_BINS > waveform = mappedData->GetTrace ();
249- int maxPos = FindMaxPosition (waveform);
250-
251- // Calculate rms before filtering
252-
253- double rmsRaw = ComputeBaselineMean (waveform, maxPos, 0 );
254- double meanRaw = ComputeBaselineMean (waveform, maxPos, 1 ) + mappedData->GetBaseline ();
255-
256- // create baseline-subtracted copy and write it to file (temporary tool)
257- std::array<double , ACTAF_BINS > waveform_bs = waveform;
258- double baseline = meanRaw;
259- for (auto & x : waveform_bs)
260- x -= baseline;
261163
262164 // Apply the SG filter to the waveform (on baseline-subtracted data)
263165 if (fApplySGFilter )
264- ApplySGFilter (waveform, fSgCoeffs );
265-
266- // Use baseline-subtracted waveform for further processing
267- // waveform = waveform_bs;
268-
269- auto integral = IntegratePulse (waveform, mappedData->GetMaxpos ());
270- auto energy = integral * fEGain [pad - 1 ];
271- auto energyMaxAmpl = mappedData->GetMaxampl () * fEGain [pad - 1 ];
166+ R3BActafUtils::ApplySGFilter (waveform, fSgCoeffs );
272167
273168 auto drift = mappedData->GetLeadingEdgeTime () * fConversionCh2ns ; // in ns
274169 auto zpos = drift * fVelocity ; // in cm
275- auto syntime = drift - synTagTime; // in ns
170+ auto syntime = drift - synTagTime * fConversionCh2ns ; // in ns
276171
277- // Calculation of RMS from the waveform after filtering
278- maxPos = FindMaxPosition (waveform);
279- double rms = ComputeBaselineMean (waveform, maxPos, 0 );
280- double mean = ComputeBaselineMean (waveform, maxPos, 1 ) + mappedData->GetBaseline ();
172+ auto maxPos = R3BActafUtils::FindMaxPosition (waveform);
173+ double rms = R3BActafUtils::ComputeBaselineMean (waveform, maxPos, 0 );
174+ double mean = R3BActafUtils::ComputeBaselineMean (waveform, maxPos, 1 );
281175
282176 // MAW parameter value
283177 double maw = mappedData->GetMaw ();
284178
179+ for (auto & x : waveform)
180+ x -= mean; // baseline
181+
182+ auto integral = R3BActafUtils::IntegratePulse (waveform, maxPos);
183+ auto energy = integral * fEGain [pad - 1 ];
184+ auto MaxAmpl = R3BActafUtils::FindMaxAmplitude (waveform);
185+ auto energyMaxAmpl = MaxAmpl * fEGain [pad - 1 ];
186+
285187 if (energyMaxAmpl >= fEThr [pad - 1 ] && energy < fMaxE )
286- AddCalData (pad, energy, energyMaxAmpl, drift, zpos, syntime, waveform, rmsRaw, rms, meanRaw, mean, maw);
188+ {
189+ AddCalData (pad,
190+ energy,
191+ energyMaxAmpl,
192+ drift,
193+ zpos,
194+ syntime,
195+ waveform,
196+ mappedData->GetRms (),
197+ rms,
198+ mappedData->GetBaseline (),
199+ mean,
200+ maw);
201+ }
287202 }
288203 return ;
289204}
0 commit comments