ExpTailRemover Class Reference

#include <ExpTailRemover.h>

Inheritance diagram for ExpTailRemover:

Public Types
using	Index = unsigned int
using	Vector = std::vector< double >
using	Name = std::string
using	NameVector = std::vector< Name >

Public Member Functions
	ExpTailRemover (fhicl::ParameterSet const &ps)
	~ExpTailRemover () override=default
DataMap	update (AdcChannelData &acd) const override

Private Member Functions
void	getSignal (const Vector &qdats, double ped, double csi, Vector &qsigs) const
void	estimatepars (const Vector &qdats, double &ped, double &csi, Index &nsamFit) const
void	wlinfit (const Vector &x, const Vector &y, const Vector &e, double &slope, double &intercept) const
Private Member Functions inherited from TpcDataTool
virtual DataMap	updateTpcData (TpcData &) const
virtual DataMap	viewTpcData (const TpcData &) const
virtual int	forwardTpcData () const
Private Member Functions inherited from AdcChannelTool
virtual	~AdcChannelTool ()=default
virtual DataMap	view (const AdcChannelData &acd) const
virtual DataMap	updateMap (AdcChannelDataMap &acds) const
virtual DataMap	viewMap (const AdcChannelDataMap &acds) const
virtual bool	updateWithView () const
virtual bool	viewWithUpdate () const
virtual DataMap	beginEvent (const DuneEventInfo &) const
virtual DataMap	endEvent (const DuneEventInfo &) const
virtual DataMap	close (const DataMap *dmin=nullptr)

Private Attributes
int	m_LogLevel
Index	m_SignalFlag
Index	m_SignalIterationLimit
Name	m_SignalTool
double	m_DecayTime
NameVector	m_IncludeChannelRanges
NameVector	m_ExcludeChannelRanges
bool	m_useChannelRanges
std::vector< bool >	m_checkChannels
AdcChannelTool *	m_pSignalTool

Additional Inherited Members
Private Types inherited from AdcChannelTool
using	Index = unsigned int
Static Private Member Functions inherited from AdcChannelTool
static int	interfaceNotImplemented ()

Detailed Description

Definition at line 72 of file ExpTailRemover.h.

Member Typedef Documentation

using ExpTailRemover::Index = unsigned int

Definition at line 76 of file ExpTailRemover.h.

using ExpTailRemover::Name = std::string

Definition at line 78 of file ExpTailRemover.h.

using ExpTailRemover::NameVector = std::vector<Name>

Definition at line 79 of file ExpTailRemover.h.

using ExpTailRemover::Vector = std::vector<double>

Definition at line 77 of file ExpTailRemover.h.

Constructor & Destructor Documentation

ExpTailRemover::ExpTailRemover

(

fhicl::ParameterSet const &

ps

)

Definition at line 26 of file ExpTailRemover_tool.cc.

: m_LogLevel(ps.get<int>("LogLevel")),
  m_SignalFlag(ps.get<Index>("SignalFlag")),
  m_SignalIterationLimit(ps.get<Index>("SignalIterationLimit")),
  m_SignalTool(ps.get<string>("SignalTool")),
  m_DecayTime(ps.get<double>("DecayTime")) ,
  m_IncludeChannelRanges(ps.get<NameVector>("IncludeChannelRanges")),
  m_ExcludeChannelRanges(ps.get<NameVector>("ExcludeChannelRanges")),
  m_useChannelRanges(false),
  m_pSignalTool(nullptr) {
  const string myname = "ExpTailRemover::ctor: ";
  if ( m_SignalFlag > 3 ) {
    Index signalFlag = 2;
    cout << myname << "WARNING: Invalid SignalFlag value " << m_SignalFlag
         << " reset to " << signalFlag << "." << endl;
    m_SignalFlag = signalFlag;
  }
  if ( m_SignalFlag >= 2  ) {
    DuneToolManager* ptm = DuneToolManager::instance();
    if ( ptm == nullptr ) {
      cout << myname << "ERROR: Unable to retrieve tool manager." << endl;
    } else {
      m_pSignalTool = ptm->getShared<AdcChannelTool>(m_SignalTool);
      if ( m_pSignalTool == nullptr ) {
        cout << myname << "ERROR: Signal finding tool not found: " << m_SignalTool << endl;
      }
    }
  }
  Index nchaCheck = 0;
  if ( m_IncludeChannelRanges.size() ) {
    DuneToolManager* ptm = DuneToolManager::instance();
    const IndexRangeTool* pcrt = ptm->getShared<IndexRangeTool>("channelRanges");
    if ( pcrt == nullptr ) {
      cout << myname << "ERROR: IndexRangeTool not found: channelRanges" << endl;
    } else {
      for ( Name crn : m_IncludeChannelRanges ) {
        IndexRange ran = pcrt->get(crn);
        if ( ran.isValid() ) {
          if ( ran.end > m_checkChannels.size() ) m_checkChannels.resize(ran.end, false);
          for ( Index icha=ran.begin; icha<ran.end; ++icha ) m_checkChannels[icha] = true;
        } else {
          cout << myname << "WARNING: Ignoring invalid include channel range " << crn << endl;
        }
      }
      if ( m_ExcludeChannelRanges.size() ) {
        for ( Name crn : m_ExcludeChannelRanges ) {
          if ( crn == "all" ) {
            m_checkChannels.clear();
            break;
          }
          IndexRange ran = pcrt->get(crn);
          if ( ran.isValid() ) {
            Index end = ran.end < m_checkChannels.size() ? ran.end : m_checkChannels.size();
            for ( Index icha=ran.begin; icha<end; ++icha ) m_checkChannels[icha] = false;
          } else {
            cout << myname << "WARNING: Ignoring invalid exclude channel range " << crn << endl;
          }
        }
      }
      m_useChannelRanges = true;
      for ( bool keep : m_checkChannels ) if ( keep ) ++nchaCheck;
    }
  }
  if ( m_LogLevel >= 1 ) {
    cout << myname << "Parameters:" << endl;
    cout << myname << "              LogLevel: " << m_LogLevel << endl;
    cout << myname << "            SignalFlag: " << m_SignalFlag << endl;
    cout << myname << "  SignalIterationLimit: " << m_SignalIterationLimit << endl;
    cout << myname << "            SignalTool: " << m_SignalTool << endl;
    cout << myname << "             DecayTime: " << m_DecayTime << endl;
    cout << myname << "  IncludeChannelRanges: [";
    bool first = true;
    for ( Name crn : m_IncludeChannelRanges ) {
      if ( first ) first = false;
      else cout << ", ";
      cout  << crn;
    }
    cout << "]" << endl;
    cout << myname << "  ExcludeChannelRanges: [";
    first = true;
    for ( Name crn : m_ExcludeChannelRanges ) {
      if ( first ) first = false;
      else cout << ", ";
      cout  << crn;
    }
    cout << "]" << endl;
    if ( m_useChannelRanges ) {
      cout << myname << "Channel checking enabled for " << nchaCheck << " channel"
           << ( nchaCheck == 1 ? "" : "s") << "." << endl;
    } else {
      cout << myname << "Channel checking disabled." << endl;
    }
  }
}

ExpTailRemover::~ExpTailRemover ( )

overridedefault

Member Function Documentation

void ExpTailRemover::estimatepars ( const Vector &  qdats, double &  ped, double &  csi, Index &  nsamFit  ) const

private

void ExpTailRemover::getSignal ( const Vector &  qdats, double  ped, double  csi, Vector &  qsigs  ) const

private

DataMap ExpTailRemover::update ( AdcChannelData &  acd ) const

overridevirtual

Reimplemented from AdcChannelTool.

Definition at line 123 of file ExpTailRemover_tool.cc.

                                                        {
  const string myname = "ExpTailRemover::update: ";
  DataMap ret;

  // Save the data before tail removal.
  AdcSignalVector samples = acd.samples;
  Index nsam = samples.size();

  // Check the channel.
  if ( m_useChannelRanges ) {
    if ( acd.channel() >= m_checkChannels.size() || ! m_checkChannels[acd.channel()] ) {
      if ( m_LogLevel >= 2 ) cout << myname << "Skipping channel " << acd.channel() << endl;
      return ret;
    }
  }

  // Check input data size.
  if ( nsam < 10 ) {
    cout << myname << "WARNING: Data for channel " << acd.channel() << " has "
         << ( nsam==0 ? "no" : "too few" ) << " ticks." << endl;
    return ret.setStatus(1);
  }

  if ( m_LogLevel >= 2 ) cout << myname << "Correcting run " << acd.run() << " event " << acd.event()
                              << " channel " << acd.channel() << endl;

  // Build the initial signal selection.
  bool checkSignal = true;   // Whether to use only non-signal in fit.
  bool findSignal = false;   // Whether to find signals each iteration.
  if ( m_SignalFlag == 0 ) {
    checkSignal = false;
  } else if ( m_SignalFlag == 1  ) {
    if ( acd.signal.size() < nsam ) {
      cout << myname << "WARNING: Data is missing signal flags--padding from " << acd.signal.size()
           << " to " << nsam << " samples." << endl;
    }
  } else if ( m_SignalFlag >= 2  ) {
    if ( m_pSignalTool == nullptr ) {
      cout << myname << "WARNING: Signal-finding tool is missing. Using all signals." << endl;
      checkSignal = false;
    } else {
      findSignal = true;
    }
  }

  Index niter = 0;  // # fit iterations
  float ped = 0.0;  // Fitted pedestal
  float tau = 0.0;  // Fitted tail0
  Index nsamKeep = 0;
  Index maxiter = findSignal ? m_SignalIterationLimit : 1;
  double noise = 0.0;
  while ( niter < maxiter ) {
    // Do signal finding.
    if ( findSignal ) {
      AdcFilterVector signalLast = acd.signal;
      DataMap fret = m_pSignalTool->update(acd);
      if ( fret ) {
        cout << myname << "WARNING: Signal-finding failed for event " << acd.event()
             << " channel " << acd.channel() << endl;
        break;
      }
      if ( acd.signal == signalLast ) {
        if ( m_LogLevel >=3 ) cout << myname << "Signal is unchanged. Exiting loop." << endl;
        break;
      }
    }
    //
    // Evaluate the signal coefficients.
    SampleTailer sta(m_DecayTime);
    sta.setTail0(1.0);
    sta.setDataZero(nsam);
    DoubleVector ctau(nsam);
    copy(sta.signal().begin(), sta.signal().end(), ctau.begin());
    sta.setTail0(0.0);
    sta.setPedestal(1.0);
    sta.setDataZero(nsam);
    DoubleVector cped(nsam);
    copy(sta.signal().begin(), sta.signal().end(), cped.begin());
    sta.setPedestal(0.0);
    sta.setData(samples);
    DoubleVector cdat(nsam);
    copy(sta.signal().begin(), sta.signal().end(), cdat.begin());
    // Evaluate the K-params.
    double kdd = 0.0;
    double kdt = 0.0;
    double kdp = 0.0;
    double ktt = 0.0;
    double ktp = 0.0;
    double kpp = 0.0;
    nsamKeep = 0;
    for ( Index  isam=0; isam<nsam; ++isam ) {
      if ( checkSignal &&  isam < acd.signal.size() && acd.signal[isam] ) continue;
      double cd = cdat[isam];
      double ct = ctau[isam];
      double cp = cped[isam];
      kdd += cd*cd;
      kdt += cd*ct;
      kdp += cd*cp;
      ktt += ct*ct;
      ktp += ct*cp;
      kpp += cp*cp;
      ++nsamKeep;
    }
    // Invert matrix and solve for (ped, tau).
    double den = ktt*kpp - ktp*ktp;
    if ( den == 0.0 ) {
      if ( acd.channelStatus() == 0 || m_LogLevel >= 2 ) {
        cout << myname << "WARNING: Unable to invert K-matrix with "
             << nsamKeep << " of " << nsam << " samples--stopping iteration for channel "
             << acd.channel() << " with status " << acd.channelStatus() << "." << endl;
      }
      break;
    }
    double deninv = 1.0/den;
    tau = deninv*(kdp*ktp-kdt*kpp);
    ped = deninv*(kdt*ktp-kdp*ktt);
    double chsq = kdd + ktt*tau*tau + kpp*ped*ped + 2.0*(kdt*tau + kdp*ped + ktp*tau*ped);
    noise = 0.0;
    if ( nsamKeep > 2 && chsq > 0.0 ) {
      noise =  sqrt(chsq/(nsamKeep-2));
    }
    if ( m_LogLevel >= 3 ) cout << myname << "Iteration " << niter << ": ped, tau, noise: "
                                << ped << ", " << tau << ", " << noise
                                << " (" << nsamKeep << " samples)." << endl;
    // Use the parameters to remove tail from the original data and store
    // these tail-subtracted samples in the channel data.
    sta.setTail0(tau);
    sta.setPedestal(ped);
    sta.setData(samples);
    acd.samples = sta.signal();
    // Update iteration count.
    ++niter;
  }

  // Log result of interation.
  if ( m_LogLevel >= 2 ) {
    cout << myname << "Iteration count: " << niter << endl;
    cout << myname << "Final ped, tau: " << ped << ", " << tau << endl;
  }

  // Use the fitted
  acd.metadata["uscPedestal"] = ped;
  acd.metadata["uscTail"] = tau;
  acd.metadata["uscNoise"] = noise;
  ret.setFloat("uscPedestal", ped);
  ret.setFloat("uscTail", tau);
  ret.setFloat("uscNoise", noise);
  ret.setInt("uscNsamFit", nsamKeep);
  ret.setInt("uscNiteration", niter);

  return ret;
}

void ExpTailRemover::wlinfit ( const Vector &  x, const Vector &  y, const Vector &  e, double &  slope, double &  intercept  ) const

private

Member Data Documentation

std::vector<bool> ExpTailRemover::m_checkChannels

private

Definition at line 100 of file ExpTailRemover.h.

double ExpTailRemover::m_DecayTime

private

Definition at line 94 of file ExpTailRemover.h.

NameVector ExpTailRemover::m_ExcludeChannelRanges

private

Definition at line 96 of file ExpTailRemover.h.

NameVector ExpTailRemover::m_IncludeChannelRanges

private

Definition at line 95 of file ExpTailRemover.h.

int ExpTailRemover::m_LogLevel

private

Definition at line 90 of file ExpTailRemover.h.

AdcChannelTool* ExpTailRemover::m_pSignalTool

private

Definition at line 101 of file ExpTailRemover.h.

Index ExpTailRemover::m_SignalFlag

private

Definition at line 91 of file ExpTailRemover.h.

Index ExpTailRemover::m_SignalIterationLimit

private

Definition at line 92 of file ExpTailRemover.h.

Name ExpTailRemover::m_SignalTool

private

Definition at line 93 of file ExpTailRemover.h.

bool ExpTailRemover::m_useChannelRanges

private

Definition at line 99 of file ExpTailRemover.h.

---

The documentation for this class was generated from the following files:

• dunedataprep/dunedataprep/DataPrep/Tool/ExpTailRemover.h
• dunedataprep/dunedataprep/DataPrep/Tool/ExpTailRemover_tool.cc