WireCell::SigProc::L1SPFilter Class Reference

#include <L1SPFilter.h>

Inheritance diagram for WireCell::SigProc::L1SPFilter:

Public Member Functions
	L1SPFilter (double gain=14.0 *units::mV/units::fC, double shaping=2.2 *units::microsecond, double postgain=1.2, double ADC_mV=4096/(2000.*units::mV), double fine_time_offset=0.0 *units::microsecond, double coarse_time_offset=-8.0 *units::microsecond)
virtual	~L1SPFilter ()
virtual bool	operator() (const input_pointer &in, output_pointer &out)
	IFrameFilter interface. More...
virtual void	configure (const WireCell::Configuration &config)
	IConfigurable interface. More...
virtual WireCell::Configuration	default_configuration () const
	Optional, override to return a hard-coded default configuration. More...
void	init_resp ()
int	L1_fit (std::shared_ptr< WireCell::SimpleTrace > &newtrace, std::shared_ptr< const WireCell::ITrace > &adctrace, int start_tick, int end_tick, bool flag_shorted=false)
Public Member Functions inherited from WireCell::IFrameFilter
virtual	~IFrameFilter ()
virtual std::string	signature ()
	Set the signature for all subclasses. More...
Public Member Functions inherited from WireCell::IFunctionNode< IFrame, IFrame >
virtual	~IFunctionNode ()
virtual bool	operator() (const boost::any &anyin, boost::any &anyout)
	The calling signature: More...
virtual std::vector< std::string >	input_types ()
virtual std::vector< std::string >	output_types ()
Public Member Functions inherited from WireCell::IFunctionNodeBase
virtual	~IFunctionNodeBase ()
virtual NodeCategory	category ()
	Return the behavior category type. More...
virtual int	concurrency ()
	By default assume all subclasses are stateless. More...
Public Member Functions inherited from WireCell::INode
virtual	~INode ()
virtual void	reset ()
Public Member Functions inherited from WireCell::IComponent< INode >
virtual	~IComponent ()
Public Member Functions inherited from WireCell::Interface
virtual	~Interface ()
Public Member Functions inherited from WireCell::IConfigurable
virtual	~IConfigurable ()
Public Member Functions inherited from WireCell::IComponent< IConfigurable >
virtual	~IComponent ()

Private Attributes
Configuration	m_cfg
double	m_gain
double	m_shaping
double	m_postgain
double	m_ADC_mV
double	m_fine_time_offset
double	m_coarse_time_offset
double	m_period
linterp< double > *	lin_V
linterp< double > *	lin_W

Additional Inherited Members
Public Types inherited from WireCell::IFrameFilter
typedef std::shared_ptr< IFrameFilter >	pointer
Public Types inherited from WireCell::IFunctionNode< IFrame, IFrame >
typedef IFrame	input_type
typedef IFrame	output_type
typedef std::shared_ptr< const IFrame >	input_pointer
typedef std::shared_ptr< const IFrame >	output_pointer
typedef IFunctionNode< IFrame, IFrame >	signature_type
Public Types inherited from WireCell::IFunctionNodeBase
typedef std::shared_ptr< IFunctionNodeBase >	pointer
Public Types inherited from WireCell::INode
enum	NodeCategory {   unknown, sourceNode, sinkNode, functionNode,   queuedoutNode, joinNode, splitNode, faninNode,   fanoutNode, multioutNode, hydraNode }
Public Types inherited from WireCell::IComponent< INode >
typedef std::shared_ptr< INode >	pointer
	Access subclass facet by pointer. More...
typedef std::vector< pointer >	vector
	Vector of shared pointers. More...
Public Types inherited from WireCell::Interface
typedef std::shared_ptr< Interface >	pointer
Public Types inherited from WireCell::IComponent< IConfigurable >
typedef std::shared_ptr< IConfigurable >	pointer
	Access subclass facet by pointer. More...
typedef std::vector< pointer >	vector
	Vector of shared pointers. More...

Detailed Description

Definition at line 18 of file L1SPFilter.h.

Constructor & Destructor Documentation

L1SPFilter::L1SPFilter	(	double	gain = 14.0 * units::mV/units::fC,
		double	shaping = 2.2 * units::microsecond,
		double	postgain = 1.2,
		double	ADC_mV = 4096/(2000.*units::mV),
		double	fine_time_offset = 0.0 * units::microsecond,
		double	coarse_time_offset = -8.0 * units::microsecond
	)

Definition at line 25 of file L1SPFilter.cxx.

  : m_gain(gain)
  , m_shaping(shaping)
  , m_postgain(postgain)
  , m_ADC_mV(ADC_mV)
  , m_fine_time_offset(fine_time_offset)
  , m_coarse_time_offset(coarse_time_offset)
  , lin_V(0)
  , lin_W(0)
{
}

L1SPFilter::~L1SPFilter ( )

virtual

Definition at line 42 of file L1SPFilter.cxx.

{
  delete lin_V;
  delete lin_W;
}

Member Function Documentation

void L1SPFilter::configure ( const WireCell::Configuration &  config )

virtual

IConfigurable interface.

Implements WireCell::IConfigurable.

Definition at line 165 of file L1SPFilter.cxx.

{
    m_cfg = cfg;

    m_gain = get(cfg,"gain",m_gain);
    m_shaping = get(cfg,"shaping",m_shaping);
    m_postgain = get(cfg,"postgain", m_postgain);
    m_ADC_mV = get(cfg,"ADC_mV", m_ADC_mV);

    m_fine_time_offset = get(cfg,"fine_time_offset", m_fine_time_offset);
    m_coarse_time_offset = get(cfg,"coarse_time_offset", m_coarse_time_offset);
}

WireCell::Configuration L1SPFilter::default_configuration ( ) const

virtual

Optional, override to return a hard-coded default configuration.

Name of component providing field responses

An array holding a waveform to use as the "smearing" filter.

The tag identifying traces which represent "raw" (not deconvolved) ADC values.

The tag identifying traces which represent "signal" processed (deconvolved) waveforms.

The tag to place on the output waveforms

Reimplemented from WireCell::IConfigurable.

Definition at line 97 of file L1SPFilter.cxx.

{
    Configuration cfg;

    /// Name of component providing field responses
    cfg["fields"] = "FieldResponse";

    /// An array holding a waveform to use as the "smearing" filter.
    cfg["filter"] = Json::arrayValue;

    /// The tag identifying traces which represent "raw" (not
    /// deconvolved) ADC values.
    cfg["adctag"] = "raw";

    /// The tag identifying traces which represent "signal" processed
    /// (deconvolved) waveforms.
    cfg["sigtag"] = "gauss";

    /// The tag to place on the output waveforms
    cfg["outtag"] = "l1sp";

    // 4 sigma for raw waveform ROI identification
    cfg["raw_ROI_th_nsigma"] = 4;
    // 10 ADC for upper limit on ADC ... 
    cfg["raw_ROI_th_adclimit"] = 10;
    // global offset 
    cfg["overall_time_offset"] = 0;
    // need 3 us offset for collection plane relative to the induction plane ...
    cfg["collect_time_offset"] = 3.0;

    // ROI padding ticks ...
    cfg["roi_pad"] = 3;
    cfg["raw_pad"] = 15;

    // L1 fit parameters ...
    cfg["adc_l1_threshold"] = 6;
    cfg["adc_sum_threshold"] = 160;
    cfg["adc_sum_rescaling"] = 90.;
    cfg["adc_sum_rescaling_limit"] = 50.;
    cfg["adc_ratio_threshold"] = 0.2;
    
    cfg["l1_seg_length"] = 120;
    cfg["l1_scaling_factor"] = 500;
    cfg["l1_lambda"] = 5;
    cfg["l1_epsilon"] = 0.05;
    cfg["l1_niteration"] = 100000;
    cfg["l1_decon_limit"] = 100; // 100 electrons

    cfg["l1_resp_scale"] = 0.5;
    cfg["l1_col_scale"] = 1.15;
    cfg["l1_ind_scale"] = 0.5;

    cfg["peak_threshold"] = 1000;
    cfg["mean_threshold"] = 500;
    
    
    cfg["gain"] = m_gain;
    cfg["shaping"] = m_shaping;
    cfg["postgain"] = m_postgain;
    cfg["ADC_mV"] = m_ADC_mV;

    cfg["fine_time_offset"] = m_fine_time_offset;
    cfg["coarse_time_offset"] = m_coarse_time_offset;

    
    return cfg;
}

void L1SPFilter::init_resp

(

)

Definition at line 48 of file L1SPFilter.cxx.

                          {
  if (lin_V==0 && lin_W==0){
   // get field response ... 
    auto ifr = Factory::find_tn<IFieldResponse>(get<std::string>(m_cfg, "fields", "FieldResponse"));
    Response::Schema::FieldResponse fr = ifr->field_response();
    // Make a new data set which is the average FR, make an average for V and Y planes ...
    Response::Schema::FieldResponse fravg = Response::average_1D(fr);

    //get electronics response
    WireCell::Waveform::compseq_t elec;
    WireCell::Binning tbins(Response::as_array(fravg.planes[0]).cols(), 0, Response::as_array(fravg.planes[0]).cols() * fravg.period);
    Response::ColdElec ce(m_gain, m_shaping);
    auto ewave = ce.generate(tbins);
    Waveform::scale(ewave, m_postgain * m_ADC_mV * (-1)); // ADC to electron ... 
    elec = Waveform::dft(ewave);

    std::complex<float> fine_period(fravg.period,0);
    
    // do a convolution here ...
    WireCell::Waveform::realseq_t resp_V = fravg.planes[1].paths[0].current ;
    WireCell::Waveform::realseq_t resp_W = fravg.planes[2].paths[0].current ; 
    
    auto spectrum_V = WireCell::Waveform::dft(resp_V);
    auto spectrum_W = WireCell::Waveform::dft(resp_W);

    WireCell::Waveform::scale(spectrum_V,elec);
    WireCell::Waveform::scale(spectrum_W,elec);
    
    WireCell::Waveform::scale(spectrum_V,fine_period);
    WireCell::Waveform::scale(spectrum_W,fine_period);
    
    // Now this response is ADC for 1 electron .
    resp_V = WireCell::Waveform::idft(spectrum_V);
    resp_W = WireCell::Waveform::idft(spectrum_W);
    
    // convolute with V and Y average responses ... 
    double intrinsic_time_offset = fravg.origin/fravg.speed;
    //std::cout << intrinsic_time_offset << " " << m_fine_time_offset << " " << m_coarse_time_offset << " " << m_gain << " " << 14.0 * units::mV/units::fC << " " << m_shaping << " " << fravg.period << std::endl;

    double x0 = (- intrinsic_time_offset - m_coarse_time_offset + m_fine_time_offset);
    double xstep = fravg.period;

    // boost::math::cubic_b_spline<double> spline_v(resp_V.begin(), resp_V.end(), x0, xstep);
    // boost::math::cubic_b_spline<double> spline_w(resp_W.begin(), resp_W.end(), x0, xstep);
    lin_V = new linterp<double>(resp_V.begin(), resp_V.end(), x0, xstep);
    lin_W = new linterp<double>(resp_W.begin(), resp_W.end(), x0, xstep);
  }
}

int L1SPFilter::L1_fit	(	std::shared_ptr< WireCell::SimpleTrace > &	newtrace,
		std::shared_ptr< const WireCell::ITrace > &	adctrace,
		int	start_tick,
		int	end_tick,
		bool	flag_shorted = false
	)

Definition at line 494 of file L1SPFilter.cxx.

                                                                                                                                                                  {
  
  double overall_time_offset = get(m_cfg,"overall_time_offset",overall_time_offset) * units::us;
  double collect_time_offset = get(m_cfg,"collect_time_offset",collect_time_offset) * units::us;
  double adc_l1_threshold = get(m_cfg,"adc_l1_threshold",adc_l1_threshold);
  double adc_sum_threshold= get(m_cfg,"adc_sum_threshold",adc_sum_threshold);
  double adc_sum_rescaling= get(m_cfg,"adc_sum_rescaling",adc_sum_rescaling);
  double adc_sum_rescaling_limit= get(m_cfg,"adc_sum_rescaling_limit",adc_sum_rescaling_limit);
  double adc_ratio_threshold = get(m_cfg,"adc_ratio_threshold",adc_ratio_threshold);
  
  double l1_seg_length= get(m_cfg,"l1_seg_length",l1_seg_length);
  double l1_scaling_factor= get(m_cfg,"l1_scaling_factor",l1_scaling_factor);
  double l1_lambda= get(m_cfg,"l1_lambda",l1_lambda);
  double l1_epsilon= get(m_cfg,"l1_epsilon",l1_epsilon);
  double l1_niteration= get(m_cfg,"l1_niteration",l1_niteration);
  double l1_decon_limit= get(m_cfg,"l1_decon_limit",l1_decon_limit);
  
  double l1_resp_scale = get(m_cfg,"l1_resp_scale",l1_resp_scale);
  double l1_col_scale = get(m_cfg,"l1_col_scale",l1_col_scale);
  double l1_ind_scale = get(m_cfg,"l1_ind_scale",l1_ind_scale);

  double mean_threshold = get(m_cfg,"mean_threshold",mean_threshold);
  double peak_threshold = get(m_cfg,"peak_threshold",peak_threshold);
  
  
  std::vector<double> smearing_vec = get< std::vector<double> >(m_cfg, "filter");
  
  // algorithm 
  const int nbin_fit = end_tick-start_tick;

  // fill the data ... 
  VectorXd init_W = VectorXd::Zero(nbin_fit);
  VectorXd init_beta = VectorXd::Zero(nbin_fit);
  VectorXd final_beta = VectorXd::Zero(nbin_fit*2);

  double temp_sum = 0;
  double temp1_sum = 0;
  double temp2_sum = 0;
  double max_val = -1;
  double min_val = 1;
  for (int i=0; i!= nbin_fit; i++){
    init_W(i) =  adctrace->charge().at(i+start_tick-newtrace->tbin()) ;
    init_beta(i) = newtrace->charge().at(i+start_tick-newtrace->tbin()) ;

    if (max_val < init_W(i)) max_val = init_W(i);
    if (min_val > init_W(i)) min_val = init_W(i);
    
    if (fabs(init_W(i))>adc_l1_threshold) {
      temp_sum += init_W(i);
      temp1_sum += fabs(init_W(i));
      temp2_sum += fabs(init_beta(i));
    }
  }

  

  int flag_l1 = 0; // do nothing
  // 1 do L1 
  if (temp_sum/(temp1_sum*adc_sum_rescaling*1.0/nbin_fit)>adc_ratio_threshold && temp1_sum>adc_sum_threshold){
    flag_l1 = 1; // do L1 ...
  }else if (temp1_sum*adc_sum_rescaling*1.0/nbin_fit < adc_sum_rescaling_limit){
    flag_l1 = 2; //remove signal ... 
  }else if (temp2_sum > 30 * nbin_fit && temp1_sum < 2.0*nbin_fit && max_val - min_val < 22){
    flag_l1 = 2;
  }

 

  // if (adctrace->channel() == 4079){
  // std::cout << adctrace->channel() << " " << nbin_fit << " " << start_tick << " " << end_tick << " " << temp_sum << " " << temp1_sum << " " << temp2_sum << " " << max_val << " " << min_val << " " << flag_l1 << std::endl;
  // }
  
  // std::cout << temp_sum << " " << temp1_sum << " " << temp_sum/(temp1_sum*adc_sum_rescaling*1.0/nbin_fit) << " " << adc_ratio_threshold << " " << temp1_sum*adc_sum_rescaling*1.0/nbin_fit << " " << flag_l1 << std::endl;
  
  if (flag_l1==1){
    // do L1 fit ...
    int n_section = std::round(nbin_fit/l1_seg_length*1.0);
    if (n_section ==0) n_section =1;
    std::vector<int> boundaries;
    for (int i=0;i!=n_section;i++){
      boundaries.push_back(int(i * nbin_fit /n_section));
    }
    boundaries.push_back(nbin_fit);

    for (int i=0;i!=n_section;i++){
      int temp_nbin_fit = boundaries.at(i+1)-boundaries.at(i);
      VectorXd W = VectorXd::Zero(temp_nbin_fit);
      for (int j=0;j!=temp_nbin_fit;j++){
        W(j) = init_W(boundaries.at(i)+j);
      }

      //for matrix G
      MatrixXd G = MatrixXd::Zero(temp_nbin_fit, temp_nbin_fit*2);

      for (int i=0;i!=temp_nbin_fit;i++){ 
        // X 
        double t1 = i/2.*units::us; // us, measured time  
        for (int j=0;j!=temp_nbin_fit;j++){
          // Y ... 
          double t2 = j/2.*units::us; // us, real signal time
          double delta_t = t1 - t2;
          if (delta_t >-15*units::us - overall_time_offset && delta_t < 10*units::us - overall_time_offset){
            G(i,j) = (*lin_W)(delta_t+overall_time_offset+collect_time_offset) * l1_scaling_factor * l1_resp_scale;
            G(i,temp_nbin_fit+j) = (*lin_V)(delta_t+overall_time_offset) * l1_scaling_factor * l1_resp_scale; 
          }
        }
      }

      double lambda = l1_lambda;//1/2.;
      WireCell::LassoModel m2(lambda, l1_niteration, l1_epsilon);
      m2.SetData(G, W);
      m2.Fit();
      VectorXd beta = m2.Getbeta();
      for (int j=0;j!=temp_nbin_fit;j++){
        final_beta(boundaries.at(i)+j) = beta(j);
        final_beta(nbin_fit + boundaries.at(i)+j) = beta(temp_nbin_fit+j);
      }
    }

    double sum1 = 0;
    double sum2 = 0;
    for (int i=0;i!=nbin_fit;i++){
      sum1 += final_beta(i);
      sum2 += final_beta(nbin_fit+i);
    }

    if (sum1 > adc_l1_threshold){
      Waveform::realseq_t l1_signal(nbin_fit,0);
      Waveform::realseq_t l2_signal(nbin_fit,0);
      for (int j=0;j!=nbin_fit;j++){
        l1_signal.at(j) = final_beta(j) * l1_col_scale + final_beta(nbin_fit+j) * l1_ind_scale;
      }
      int mid_bin = (smearing_vec.size()-1)/2;
      //std::cout << smearing_vec.size() << " " << mid_bin << std::endl;
      for (int j=0;j!=nbin_fit;j++){
        double content = l1_signal.at(j);
        if (content>0){
          for (size_t k=0;k!=smearing_vec.size();k++){
            int bin = j+k-mid_bin;
            if (bin>=0&&bin<nbin_fit)
              l2_signal.at(bin) += content * smearing_vec.at(k);
          }
        }
      }
      
      for (int j=0;j!=nbin_fit;j++){
        if (l2_signal.at(j)<l1_decon_limit/l1_scaling_factor){ // 50 electrons
          l1_signal.at(j) = 0;
        }else{
          l1_signal.at(j) = l2_signal.at(j) * l1_scaling_factor; 
        }
      }

      {
        // go through the new data and clean the small peaks ...
        std::vector<int> nonzero_bins;
        std::vector<std::pair<int,int>> ROIs;
        for (int j=0;j!=nbin_fit;j++){
          if (l1_signal.at(j)>0)
            nonzero_bins.push_back(j);
        }

        bool flag_start = false;
        int start_bin=-1,end_bin=-1;
        for (size_t j=0;j!=nonzero_bins.size();j++){
          if (!flag_start){
            start_bin = nonzero_bins.at(j);
            end_bin = nonzero_bins.at(j);
            flag_start = true;
          }else{
            if (nonzero_bins.at(j) - end_bin == 1){
              end_bin = nonzero_bins.at(j);
            }else{
              ROIs.push_back(std::make_pair(start_bin,end_bin));
              start_bin = nonzero_bins.at(j);
              end_bin = nonzero_bins.at(j);
            }
          }
        }
        if (start_bin >=0)
          ROIs.push_back(std::make_pair(start_bin,end_bin));

        for (size_t j=0;j!=ROIs.size();j++){
          double max_val = -1;
          double mean_val1 = 0;
          double mean_val2 = 0;
          for (int k=ROIs.at(j).first; k<=ROIs.at(j).second; k++){
            if (l1_signal.at(k) > max_val) max_val = l1_signal.at(k);
            mean_val1 += l1_signal.at(k);
            mean_val2 ++;
          }
          if (mean_val2!=0)
            mean_val1 /= mean_val2;
          if (max_val < peak_threshold && mean_val1 < mean_threshold){
            for (int k=ROIs.at(j).first; k<=ROIs.at(j).second; k++){
              l1_signal.at(k) = 0;
            }
          }
          //std::cout << max_val << " " << mean_val1 << std::endl;
          
        }
        
        // std::cout << nonzero_bins.front() << " X " << nonzero_bins.back() << std::endl;
        // std::cout << ROIs.size() << std::endl;
        // for (size_t i=0;i!=ROIs.size();i++){
        //   std::cout << ROIs.at(i).first << " " << ROIs.at(i).second << std::endl;
        // }

        // finish cleaning ... 
      }

      
      for (int time_tick = start_tick; time_tick!= end_tick; time_tick++){
        newtrace->charge().at(time_tick-newtrace->tbin())=l1_signal.at(time_tick-start_tick);
      }
    }
  }else if (flag_l1==2){
    if (flag_shorted){
      for (int time_tick = start_tick; time_tick!= end_tick; time_tick++){
        // temporary hack to reset the data ... 
        newtrace->charge().at(time_tick-newtrace->tbin())=0;
      }
    }
  }

  return flag_l1;
  
}

bool L1SPFilter::operator() ( const input_pointer &  in, output_pointer &  out  )

virtual

IFrameFilter interface.

here, use the ADC and signal traces to do L1SP put result in out_traces

Implements WireCell::IFunctionNode< IFrame, IFrame >.

Definition at line 178 of file L1SPFilter.cxx.

{
    out = nullptr;
    if (!in) {
        return true;            // eos
    }

    std::string adctag = get<std::string>(m_cfg, "adctag");
    std::string sigtag = get<std::string>(m_cfg, "sigtag");
    std::string outtag = get<std::string>(m_cfg, "outtag");
    

    //    std::cout << smearing_vec.size() << std::endl;
    
    int roi_pad = 0;
    roi_pad = get(m_cfg,"roi_pad",roi_pad);
    int raw_pad = 0;
    raw_pad = get(m_cfg,"raw_pad",raw_pad);
    
    double raw_ROI_th_nsigma = get(m_cfg,"raw_ROI_th_nsigma",raw_ROI_th_nsigma);
    double raw_ROI_th_adclimit = get(m_cfg,"raw_ROI_th_adclimit",raw_ROI_th_adclimit);
    
    
    // std::cout << "Xin: " << raw_ROI_th_nsigma << " " << raw_ROI_th_adclimit << " " << overall_time_offset << " " << collect_time_offset << " " << roi_pad << " " << adc_l1_threshold << " " << adc_sum_threshold << " " << adc_sum_rescaling << " " << adc_sum_rescaling_limit << " " << l1_seg_length << " " << l1_scaling_factor << " " << l1_lambda << " " << l1_epsilon << " " << l1_niteration << " " << l1_decon_limit << " " << l1_resp_scale << " " << l1_col_scale << " " << l1_ind_scale << std::endl;
    init_resp();

    
   

    // std::cout << (*lin_V)(0*units::us) << " " << (*lin_W)(0*units::us) << std::endl;
    // std::cout << (*lin_V)(1*units::us) << " " << (*lin_W)(1*units::us) << std::endl;
    //    for (size_t i=0; i!=resp_V.size(); i++){
    // std::cout << (i*fravg.period - intrinsic_time_offset - m_coarse_time_offset + m_fine_time_offset)/units::us << " " << resp_V.at(i) << " " << resp_W.at(i) << " " << ewave.at(i) << std::endl;
    //}
    // std::complex<float> fine_period(fravg.period,0);
    // int fine_nticks = Response::as_array(fravg.planes[0]).cols();
    //    Waveform::realseq_t ftbins(fine_nticks);
    // for (int i=0;i!=fine_nticks;i++){
    //  ftbins.at(i) = i * fravg.period;
    //}


    
    auto adctraces = FrameTools::tagged_traces(in, adctag);
    auto sigtraces = FrameTools::tagged_traces(in, sigtag);

    if (adctraces.empty() or sigtraces.empty() or adctraces.size() != sigtraces.size()) {
        std::cerr << "L1SPFilter got unexpected input: "
                  << adctraces.size() << " ADC traces and "
                  << sigtraces.size() << " signal traces\n";
        THROW(RuntimeError() << errmsg{"L1SPFilter: unexpected input"});
    }

    m_period = in->tick();

    //std::cout << m_period/units::us << std::endl;
    
    /// here, use the ADC and signal traces to do L1SP
    ///  put result in out_traces
    ITrace::vector out_traces;
    
    std::map<int,std::set<int>> init_map;
    // do ROI from the decon signal
    for (auto trace : sigtraces) {
      int ch = trace->channel();
      int tbin = trace->tbin();
      auto const& charges = trace->charge();
      const int ntbins = charges.size();
      std::set<int> time_ticks;

      for (int qi = 0; qi < ntbins; qi++){
        if (charges[qi]>0){
          time_ticks.insert(tbin+qi);
        }
      }
      
      init_map[ch] = time_ticks;
      // if (time_ticks.size()>0){
      //        std::cout << ch << " " << time_ticks.size() << std::endl;
      // }
    }
    
    // do ROI from the raw signal
    int ntot_ticks=0;

    std::map<int,std::shared_ptr<const WireCell::ITrace>> adctrace_ch_map;
    
    for (auto trace : adctraces) {
      int ch = trace->channel();
      
      adctrace_ch_map[ch] = trace;

      int tbin = trace->tbin();
      auto const& charges = trace->charge();
      const int ntbins = charges.size();
      std::set<int>& time_ticks = init_map[ch];

      if (ntot_ticks < ntbins)
        ntot_ticks = ntbins;
      
      Waveform::realseq_t tmp_charge = charges;
      double mean = Waveform::percentile(tmp_charge,0.5);
      double mean_p1sig = Waveform::percentile(tmp_charge,0.5+0.34);
      double mean_n1sig = Waveform::percentile(tmp_charge,0.5-0.34);
      double cut = raw_ROI_th_nsigma * sqrt((pow(mean_p1sig-mean,2)+pow(mean_n1sig-mean,2))/2.);
      if (cut < raw_ROI_th_adclimit) cut = raw_ROI_th_adclimit;

      // if (ch==4090)
      //        std::cout << cut << " " << raw_pad << std::endl;
      
      for (int qi = 0; qi < ntbins; qi++){
        if (fabs(charges[qi])>cut){
          for (int qii = -raw_pad; qii!=raw_pad+1;qii++){
            if (tbin+qi+qii >=0 && tbin+qi+qii<ntot_ticks)
              time_ticks.insert(tbin+qi+qii);
          }
        }
      }
      // if (time_ticks.size()>0){
      //        std::cout << ch << " " << time_ticks.size() << std::endl;
      // }
    }

    
    // create ROIs ... 
    std::map<int, std::vector<std::pair<int,int>>> map_ch_rois;
    
    for (auto it = init_map.begin(); it!=init_map.end(); it++){
      int wire_index = it->first;
      std::set<int>& time_slices_set = it->second;
      if (time_slices_set.size()==0) continue;
      std::vector<int> time_slices;
      std::copy(time_slices_set.begin(), time_slices_set.end(), std::back_inserter(time_slices));
      
      std::vector<std::pair<int,int>> rois;
      std::vector<std::pair<int,int>> rois_save;
      
      rois.push_back(std::make_pair(time_slices.front(),time_slices.front()));
      for (size_t i=1; i<time_slices.size();i++){
        if (time_slices.at(i) - rois.back().second <= roi_pad*2){
          rois.back().second = time_slices.at(i);
        }else{
          rois.push_back(std::make_pair(time_slices.at(i),time_slices.at(i)));
        }
      }
      
      // extend the rois to both side according to the bin content
      for (auto it = rois.begin(); it!= rois.end();  it++){
        int start_bin = it->first;
        int end_bin = it->second;
        start_bin = start_bin - roi_pad;
        end_bin = end_bin + roi_pad;
        if (start_bin <0) start_bin = 0;
        if (end_bin>ntot_ticks-1) end_bin = ntot_ticks-1;
        it->first = start_bin;
        it->second = end_bin;
      }
      
      for (auto it = rois.begin(); it!= rois.end();  it++){
        if (rois_save.size()==0){
          rois_save.push_back(*it);
        }else if (it->first <= rois_save.back().second){
          rois_save.back().second = it->second;
        }else{
          rois_save.push_back(*it);
        }
      }

      if (rois_save.size()>0)
        map_ch_rois[wire_index] = rois_save;
      // for (auto it = rois_save.begin(); it!=rois_save.end(); it++){
      // std::cout << wire_index << " " << it->first << " " << it->second +1 << std::endl;
      // }
    }
    

    std::map<int, std::vector<int> > map_ch_flag_rois;
    // prepare for the output signal ...
    for (auto trace : sigtraces) {
      auto newtrace = std::make_shared<SimpleTrace>(trace->channel(), trace->tbin(), trace->charge());
      // How to access the sigtraces together ???
      if (map_ch_rois.find(trace->channel()) != map_ch_rois.end()){
        std::vector<std::pair<int,int>>& rois_save = map_ch_rois[trace->channel()];
        std::vector<int> flag_rois;
        map_ch_flag_rois[trace->channel()] = flag_rois;
        
        bool flag_shorted = false;
        for (size_t i1 = 0; i1!=rois_save.size(); i1++){
          //
                  
          int temp_flag = L1_fit(newtrace, adctrace_ch_map[trace->channel()], rois_save.at(i1).first, rois_save.at(i1).second+1, flag_shorted);
          map_ch_flag_rois[trace->channel()].push_back(temp_flag);
          
          for (int time_tick = rois_save.at(i1).first; time_tick!=rois_save.at(i1).second+1; time_tick++){
          // temporary hack to reset the data ...
            if (newtrace->charge().at(time_tick-trace->tbin())<0)
              newtrace->charge().at(time_tick-trace->tbin())=0;
          }
        }
        
        
        flag_shorted = false;
        int prev_time_tick = -2000;
        for (size_t i1 = 0; i1!=rois_save.size(); i1++){
          if (rois_save.at(i1).first - prev_time_tick > 20) flag_shorted = false;
           
           if (map_ch_flag_rois[trace->channel()].at(i1)==1){
             flag_shorted = true;
           }else if (map_ch_flag_rois[trace->channel()].at(i1)==2){
             if (flag_shorted){
               L1_fit(newtrace, adctrace_ch_map[trace->channel()], rois_save.at(i1).first, rois_save.at(i1).second+1, flag_shorted);
               map_ch_flag_rois[trace->channel()].at(i1) = 0;
             }
           }else if (map_ch_flag_rois[trace->channel()].at(i1)==0){
             flag_shorted = false;
           }
           prev_time_tick = rois_save.at(i1).second;
        }
        
        if (rois_save.size()>0){
          prev_time_tick = rois_save.back().second + 2000;

          for (size_t i1 = 0; i1!=rois_save.size(); i1++){
            if (prev_time_tick - rois_save.at(rois_save.size()-1-i1).second > 20) flag_shorted = false;
            
            if (map_ch_flag_rois[trace->channel()].at(i1)==1){
              flag_shorted = true;
            }else if (map_ch_flag_rois[trace->channel()].at(i1)==2){
              if (flag_shorted){
                L1_fit(newtrace, adctrace_ch_map[trace->channel()], rois_save.at(i1).first, rois_save.at(i1).second+1, flag_shorted);
                map_ch_flag_rois[trace->channel()].at(i1) = 0;
              }
            }else if (map_ch_flag_rois[trace->channel()].at(i1)==0){
              flag_shorted = false;
            }
            prev_time_tick = rois_save.at(i1).first;
          }
        }
        

        
      }

      
      // std::cout << trace->channel() << std::endl;
      out_traces.push_back(newtrace);
    }

    for (size_t i2 = 0; i2!=out_traces.size(); i2++){
      auto new_trace = std::make_shared<SimpleTrace>(out_traces.at(i2)->channel(),  out_traces.at(i2)->tbin(),  out_traces.at(i2)->charge());
      int ch = out_traces.at(i2)->channel();
      if (map_ch_rois.find(ch) != map_ch_rois.end()){
        std::vector<std::pair<int,int>>& rois_save = map_ch_rois[ch];
        std::vector<std::pair<int,int> > next_rois_save;
        if (map_ch_rois.find(ch+1)!=map_ch_rois.end())
          next_rois_save = map_ch_rois[ch+1];
        std::vector<std::pair<int,int> > prev_rois_save;
        if (map_ch_rois.find(ch-1)!=map_ch_rois.end())
          prev_rois_save = map_ch_rois[ch-1];

        for (size_t i1 = 0; i1!=rois_save.size(); i1++){
          if ( map_ch_flag_rois[ch].at(i1) == 2 ){
            bool flag_shorted = false;
            // need to check now ...
            for (size_t i3 = 0; i3 != next_rois_save.size(); i3++){
              if (map_ch_flag_rois[ch+1].at(i3)==1){
                if (rois_save.at(i1).first - 3 <= next_rois_save.at(i3).second + 3 &&
                    rois_save.at(i1).second + 3 >= next_rois_save.at(i3).first - 3){
                  flag_shorted = true;
                  break;
                }
              }
            }
            if (!flag_shorted)
              for (size_t i3 = 0; i3 != prev_rois_save.size(); i3++){
                if (map_ch_flag_rois[ch-1].at(i3)==1){
                  if (rois_save.at(i1).first - 3 <= prev_rois_save.at(i3).second + 3 &&
                      rois_save.at(i1).second + 3 >= prev_rois_save.at(i3).first - 3){
                    flag_shorted = true;
                    break;
                  }
                }
              }
            if(flag_shorted){
              for (int time_tick = rois_save.at(i1).first; time_tick!=rois_save.at(i1).second+1; time_tick++){
                new_trace->charge().at(time_tick-out_traces.at(i2)->tbin())=0;
              }
            }
          }
        }
      }
      out_traces.at(i2) = new_trace;
    }



    std::cerr << "L1SPFilter: frame: " << in->ident() << " "
              << adctag << "[" << adctraces.size() << "] + "
              << sigtag << "[" << sigtraces.size() << "] --> "
              << outtag << "[" << out_traces.size() << "]\n";


    
    // Finally, we save the traces to an output frame with tags.

    IFrame::trace_list_t tl(out_traces.size());
    std::iota(tl.begin(), tl.end(), 0);
    
    auto sf = new SimpleFrame(in->ident(), in->time(), out_traces, in->tick());
    sf->tag_traces(outtag, tl);
    out = IFrame::pointer(sf);
    
    return true;
}

Member Data Documentation

linterp<double>* WireCell::SigProc::L1SPFilter::lin_V

private

Definition at line 52 of file L1SPFilter.h.

linterp<double>* WireCell::SigProc::L1SPFilter::lin_W

private

Definition at line 53 of file L1SPFilter.h.

double WireCell::SigProc::L1SPFilter::m_ADC_mV

private

Definition at line 47 of file L1SPFilter.h.

Configuration WireCell::SigProc::L1SPFilter::m_cfg

private

Definition at line 42 of file L1SPFilter.h.

double WireCell::SigProc::L1SPFilter::m_coarse_time_offset

private

Definition at line 49 of file L1SPFilter.h.

double WireCell::SigProc::L1SPFilter::m_fine_time_offset

private

Definition at line 48 of file L1SPFilter.h.

double WireCell::SigProc::L1SPFilter::m_gain

private

Definition at line 44 of file L1SPFilter.h.

double WireCell::SigProc::L1SPFilter::m_period

private

Definition at line 50 of file L1SPFilter.h.

double WireCell::SigProc::L1SPFilter::m_postgain

private

Definition at line 46 of file L1SPFilter.h.

double WireCell::SigProc::L1SPFilter::m_shaping

private

Definition at line 45 of file L1SPFilter.h.

---

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

• wire-cell-build/sigproc/inc/WireCellSigProc/L1SPFilter.h
• wire-cell-build/sigproc/src/L1SPFilter.cxx