218 const string myname =
"ExpTailPedRemover::update: ";
224 Index nsam = samples.size();
238 Index ncof = ntai + nped;
240 cout << myname <<
"WARNING: Data for channel " << acd.
channel() <<
" has " 241 << ( nsam==0 ?
"no" :
"too few" ) <<
" ticks: " << nsam <<
" < " << ncof <<
endl;
242 return ret.setStatus(1);
245 cout << myname <<
"WARNING: Data for channel " << acd.
channel() <<
" has too many ticks:" 246 << nsam <<
" > " <<
m_MaxTick <<
". Please increase MaxTick." <<
endl;
247 return ret.setStatus(1);
250 if (
m_LogLevel >= 3 ) cout << myname <<
"Correcting run " << acd.
run() <<
" event " << acd.
event()
255 bool checkSignalDefault =
true;
256 bool findSignal =
false;
258 checkSignalDefault =
false;
260 if ( acd.
signal.size() < nsam ) {
261 cout << mychan <<
"WARNING: Data is missing signal flags--padding from " << acd.
signal.size()
262 <<
" to " << nsam <<
" samples." <<
endl;
266 cout << mychan <<
"WARNING: Signal-finding tool is missing. Using all signals." <<
endl;
267 checkSignalDefault =
false;
280 while ( niter < maxiter ) {
286 cout << myname <<
"WARNING: Signal-finding failed for event " << acd.
event()
290 if ( acd.
signal == signalLast && niter > 0 ) {
291 if (
m_LogLevel >=3 ) cout << mychan <<
"Signal is unchanged. Exiting loop." <<
endl;
296 bool checkSignal = checkSignalDefault;
299 for (
Index isam=0; isam<nsam; ++isam ) {
300 if ( isam < acd.
signal.size() && acd.
signal[isam] )
continue;
301 if ( ++nchk >= ncof )
break;
305 string sstat = cstat==0 ?
"Good" : cstat==1 ?
"Bad" : cstat==2 ?
"Noisy" :
"Unknown";
307 if ( dowarn ) cout << myname <<
"WARNING: " << sstat <<
" channel " << acd.
channel()
308 <<
": Not-signal sample count of " << nchk
309 <<
" is not sufficient to evaluate the " << ncof <<
" fit parameters.";
311 if ( dowarn ) cout <<
" Using all samples." <<
endl;
314 if ( dowarn ) cout <<
" Exiting loop." <<
endl;;
321 using DoubleVectorVector = std::vector<DoubleVector>;
322 DoubleVectorVector cofvecs(ncof,
DoubleVector(nsam, 0.0));
324 if ( !
useTail() ) sta.setBeta(1.0,
true);
325 if ( ! sta.isValid() ) {
326 cout << mychan <<
"ERROR: SampleTailer is invalid. Exiting." <<
endl;
333 sta.setPedestal(0.0);
334 sta.setDataZero(nsam);
336 copy(sta.signal().begin(), sta.signal().end(), ctau.begin());
339 for (
Index iped=0; iped<nped; ++iped ) {
342 sta.setDataZero(nsam);
344 copy(sta.signal().begin(), sta.signal().end(), cped.begin());
348 sta.setPedestal(0.0);
349 sta.setData(samples);
351 copy(sta.signal().begin(), sta.signal().end(), cdat.begin());
353 TMatrixDSym kmat(ncof);
356 for (
Index isam=0; isam<nsam; ++isam ) {
357 if ( checkSignal && isam < acd.
signal.size() && acd.
signal[isam] )
continue;
358 double cd = cdat[isam];
359 for ( icof=0; icof<ncof; ++icof ) {
360 double ci = cofvecs[icof][isam];
362 for (
Index jcof=0; jcof<=icof; ++jcof ) {
363 double cj = cofvecs[jcof][isam];
364 kmat[icof][jcof] += ci*cj;
370 for ( icof=0; icof<ncof; ++icof ) {
371 for (
Index jcof=0; jcof<icof; ++jcof ) {
372 kmat[jcof][icof] = kmat[icof][jcof];
378 if ( ! kmat.IsValid() || det == 0.0 ) {
380 cout << myname <<
"WARNING: Channel " << acd.
channel() <<
": Unable to invert K-matrix with " 381 << nsamFit <<
" of " << nsam <<
" samples--stopping iteration for channel " 383 for ( icof=0; icof<ncof; ++icof ) {
385 for (
Index jcof=0; jcof<ncof; ++jcof ) {
386 cout <<
setw(20) << kmat[icof][jcof];
395 for ( icof=0; icof<ncof; ++icof ) {
396 cofs[icof] = kvec[icof];
400 float tau =
useTail() ? cofs[0] : 0.0;
403 for (
Index iped=0; iped<nped; ++iped ) {
404 for (
Index isam=0; isam<nsam; ++isam ) {
410 sta.setPedestalVector(&peds);
411 sta.setData(samples);
415 for (
Index isam=0; isam<nsam; ++isam ) {
416 if ( checkSignal && isam < acd.
signal.size() && acd.
signal[isam] )
continue;
417 double sig = sta.signal()[isam];
420 Index ndof = nsamFit > ncof ? nsamFit - ncof : 1;
421 noise = sqrt(sumsq/ndof);
424 cout << mychan <<
"Iteration " << niter <<
": nsamfit=" << nsamFit
427 for (
float cof : cofs ) {
428 if ( first ) first =
false;
439 cout << mychan <<
"Iteration count: " << niter <<
endl;
440 cout << mychan <<
"Noise: " << noise;
441 if ( sampleUnit.size() ) cout <<
" " << sampleUnit;
442 cout <<
" from " << nsamFit <<
"/" << nsam <<
" samples" <<
endl;
443 cout << mychan <<
"Final " << m_fitNamesString <<
": {";
445 for (
float cof : cofs ) {
446 if ( first ) first =
false;
454 acd.
metadata[
"uscPedestal"] = cofs[1];
457 for (
Index icof=0; icof<cofs.size(); ++icof ) {
459 ret.setFloat(name, cofs[icof]);
461 ret.setFloat(
"uscNoise", noise);
462 ret.setInt(
"uscNsamFit", nsamFit);
463 ret.setInt(
"uscNiteration", niter);
IndexSet m_nowarnStatuses
ChannelGroupService::Name Name
std::vector< double > DoubleVector
AdcChannelTool * m_pSignalTool
Q_EXPORT QTSManip setw(int w)
Index m_SignalIterationLimit
Index channelStatus() const
std::vector< bool > m_checkChannels
std::vector< bool > AdcFilterVector
FloatVectorVector m_pedVectors
std::vector< AdcSignal > AdcSignalVector
Dft::FloatVector FloatVector
std::string to_string(ModuleType const mt)
QTextStream & endl(QTextStream &s)