@@ -90,9 +90,8 @@ void TRestRawSignalRecoverSaturationProcess::Initialize() {
90
90
fProcessAllSignals = false ;
91
91
fNBinsIfNotSaturated = 20 ;
92
92
fMinSaturationValue = 0 ;
93
- fBaseLineRange = TVector2 (-1 , -1 ); // -1 means no baseline range
94
- fFitRange = TVector2 (-1 , -1 ); // -1 means no fit range
95
-
93
+ fBaseLineRange = TVector2 (-1 , -1 ); // -1 means no baseline range
94
+ fFitRange = TVector2 (-1 , -1 ); // -1 means no fit range
96
95
}
97
96
98
97
// /////////////////////////////////////////////
@@ -132,7 +131,8 @@ TRestEvent* TRestRawSignalRecoverSaturationProcess::ProcessEvent(TRestEvent* evI
132
131
133
132
// Set the baseline range if it has been provided
134
133
if (fBaseLineRange .X () != -1 && fBaseLineRange .Y () != -1 )
135
- fAnaEvent ->SetBaseLineRange (fBaseLineRange .X (), fBaseLineRange .Y ()); // this will also calculate the baseline
134
+ fAnaEvent ->SetBaseLineRange (fBaseLineRange .X (),
135
+ fBaseLineRange .Y ()); // this will also calculate the baseline
136
136
137
137
// process each signal in the event
138
138
for (int s = 0 ; s < fAnaEvent ->GetNumberOfSignals (); s++) {
@@ -168,14 +168,15 @@ TRestEvent* TRestRawSignalRecoverSaturationProcess::ProcessEvent(TRestEvent* evI
168
168
if (fProcessAllSignals && saturatedBins.size () < fMinSaturatedBins ) {
169
169
saturatedBins.clear ();
170
170
// set saturated bins around maxPeakBin
171
- for (size_t i = maxPeakBin-fNBinsIfNotSaturated /2 ; i < maxPeakBin+fNBinsIfNotSaturated /2 && i<(size_t )signal ->GetNumberOfPoints (); i++) {
171
+ for (size_t i = maxPeakBin - fNBinsIfNotSaturated / 2 ;
172
+ i < maxPeakBin + fNBinsIfNotSaturated / 2 && i < (size_t )signal ->GetNumberOfPoints (); i++) {
172
173
saturatedBins.push_back (i);
173
174
}
174
175
}
175
176
176
177
if (!saturatedBins.empty ()) {
177
- RESTDebug << " Saturated bins:" << saturatedBins.front () << " to " << saturatedBins.back () <<
178
- " at " << maxValue << RESTendl;
178
+ RESTDebug << " Saturated bins:" << saturatedBins.front () << " to " << saturatedBins.back ()
179
+ << " at " << maxValue << RESTendl;
179
180
}
180
181
181
182
// Create TGraph with the not saturated bins for the fit
@@ -202,13 +203,13 @@ TRestEvent* TRestRawSignalRecoverSaturationProcess::ProcessEvent(TRestEvent* evI
202
203
" [0]+[1]*TMath::Exp(-3. * (x-[3])/[2]) * "
203
204
" (x-[3])/[2] * (x-[3])/[2] * (x-[3])/[2] / "
204
205
" (1+TMath::Exp(-10000*(x-[3])))" ,
205
- startFitRange, endFitRange);
206
+ startFitRange, endFitRange);
206
207
RESTDebug << " nPoints" << signal ->GetNumberOfPoints () << RESTendl;
207
208
RESTDebug << " Function created" << RESTendl;
208
209
// First estimation of the parameters
209
210
auto peakposEstimate = maxPeakBin + saturatedBins.size () / 2 ;
210
211
Double_t amplEstimate = maxValue;
211
- Double_t widthEstimate = (endFitRange - startFitRange)* 0.2 ;
212
+ Double_t widthEstimate = (endFitRange - startFitRange) * 0.2 ;
212
213
Double_t baselineEstimate = 250 ;
213
214
214
215
// Second (and better) estimation of the parameters
@@ -223,30 +224,34 @@ TRestEvent* TRestRawSignalRecoverSaturationProcess::ProcessEvent(TRestEvent* evI
223
224
// the amplitude estimate should be at least the maximum value of the signal
224
225
if (amplEstimate < maxValue) amplEstimate = maxValue;
225
226
}
226
- // signal->CalculateBaseLine(20,150);
227
+ // signal->CalculateBaseLine(20,150);
227
228
if (signal ->isBaseLineInitialized ()) {
228
229
baselineEstimate = signal ->GetBaseLine ();
229
230
}
230
231
RESTDebug << " Estimations: ampl=" << amplEstimate << " width=" << widthEstimate
231
232
<< " baseline=" << baselineEstimate << " peakpos=" << peakposEstimate << RESTendl;
232
233
// Configure the fit parameters
233
234
f->SetParNames (" Baseline" , " Amplitude" , " ShapingTime" , " PeakPosition" );
234
- // f->SetParameters(baselineEstimate, amplEstimate / 0.0498, widthEstimate, peakposEstimate - widthEstimate);
235
- // Baseline
235
+ // f->SetParameters(baselineEstimate, amplEstimate / 0.0498, widthEstimate, peakposEstimate -
236
+ // widthEstimate);
237
+ // Baseline
236
238
f->SetParameter (0 , baselineEstimate);
237
- f->SetParLimits (0 , 0 , maxValue); // baseline should be positive and less than the saturation value
238
- if (signal ->isBaseLineInitialized ()) { // fix the baseline to make it faster and more reliable
239
+ f->SetParLimits (0 , 0 , maxValue); // baseline should be positive and less than the saturation value
240
+ if (signal ->isBaseLineInitialized ()) { // fix the baseline to make it faster and more reliable
239
241
f->FixParameter (0 , baselineEstimate);
240
- }// */
242
+ } // */
241
243
// Amplitude
242
- f->SetParameter (1 , amplEstimate / 0.0498 ); // 0.0498=e^{-3}
243
- f->SetParLimits (1 , 0 , maxValue / 0.0498 * 100 ); // max allowed amplitude is 100 times the saturation value
244
+ f->SetParameter (1 , amplEstimate / 0.0498 ); // 0.0498=e^{-3}
245
+ f->SetParLimits (1 , 0 ,
246
+ maxValue / 0.0498 * 100 ); // max allowed amplitude is 100 times the saturation value
244
247
// Width or shaping time
245
248
f->SetParameter (2 , widthEstimate);
246
- f->SetParLimits (2 , 0 , signal ->GetNumberOfPoints ()); // width should be positive and less than the window
249
+ f->SetParLimits (2 , 0 ,
250
+ signal ->GetNumberOfPoints ()); // width should be positive and less than the window
247
251
// Peak position
248
252
f->SetParameter (3 , peakposEstimate - widthEstimate);
249
- f->SetParLimits (3 , 0 , signal ->GetNumberOfPoints ()); // peak position should be positive and less than the window
253
+ f->SetParLimits (
254
+ 3 , 0 , signal ->GetNumberOfPoints ()); // peak position should be positive and less than the window
250
255
251
256
std::string fitOptions = " R" ;
252
257
if (GetVerboseLevel () < TRestStringOutput::REST_Verbose_Level::REST_Debug) {
@@ -263,7 +268,7 @@ TRestEvent* TRestRawSignalRecoverSaturationProcess::ProcessEvent(TRestEvent* evI
263
268
for (size_t i = 0 ; i < (size_t )signal ->GetNumberOfPoints (); i++) {
264
269
if (std::find (saturatedBins.begin (), saturatedBins.end (), i) != saturatedBins.end ()) {
265
270
Double_t value = f->Eval (i) - maxValue;
266
- if (value > 0 || fProcessAllSignals ) {
271
+ if (value > 0 || fProcessAllSignals ) {
267
272
toAddSignal.AddPoint (value);
268
273
anyBinRecovered = true ;
269
274
addedIntegral += value;
@@ -272,7 +277,7 @@ TRestEvent* TRestRawSignalRecoverSaturationProcess::ProcessEvent(TRestEvent* evI
272
277
continue ;
273
278
}
274
279
}
275
- toAddSignal.AddPoint ((Short_t) 0 );
280
+ toAddSignal.AddPoint ((Short_t)0 );
276
281
}
277
282
278
283
if (GetVerboseLevel () >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) {
0 commit comments