caldata::CalWireSBND Class Reference

Inheritance diagram for caldata::CalWireSBND:

Public Types
using	Waveform = std::vector< float >
using	CandidateROI = std::pair< size_t, size_t >
using	CandidateROIVec = std::vector< CandidateROI >

Public Member Functions
	CalWireSBND (fhicl::ParameterSet const &pset)
virtual	~CalWireSBND ()
void	produce (art::Event &evt)
void	beginJob ()
void	endJob ()
void	reconfigure (fhicl::ParameterSet const &p)

Public Attributes
std::unique_ptr < sbnd_tool::IROIFinder >	fROITool

Private Member Functions
void	SubtractBaseline (std::vector< float > &holder)
void	SubtractBaselineAdv (std::vector< float > &holder)

Private Attributes
bool	fDoBaselineSub
	subtract baseline to restore DC component post-deconvolution More...
bool	fDoAdvBaselineSub
	use interpolation-based baseline subtraction More...
int	fBaseSampleBins
	bin grouping size in "interpolate" method More...
float	fBaseVarCut
	baseline variance cut used in "interpolate" method More...
std::string	fDigitModuleLabel
	module that made digits More...
std::string	fSpillName

Detailed Description

Definition at line 52 of file CalWireSBND_module.cc.

Member Typedef Documentation

using caldata::CalWireSBND::CandidateROI = std::pair<size_t, size_t>

Definition at line 74 of file CalWireSBND_module.cc.

using caldata::CalWireSBND::CandidateROIVec = std::vector<CandidateROI>

Definition at line 75 of file CalWireSBND_module.cc.

using caldata::CalWireSBND::Waveform = std::vector<float>

Definition at line 71 of file CalWireSBND_module.cc.

Constructor & Destructor Documentation

caldata::CalWireSBND::CalWireSBND ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 101 of file CalWireSBND_module.cc.

  : EDProducer(pset)
  {
    fSpillName="";
    this->reconfigure(pset);
    
    fROITool = art::make_tool<sbnd_tool::IROIFinder>(pset.get<fhicl::ParameterSet>("ROITool"));

    //--Hec if(fSpillName.size()<1) produces< std::vector<recob::Wire> >();
    //--Hec else produces< std::vector<recob::Wire> >(fSpillName);
    produces< std::vector<recob::Wire> >(fSpillName);
    produces<art::Assns<raw::RawDigit, recob::Wire>>(fSpillName);
    
  }

caldata::CalWireSBND::~CalWireSBND ( )

virtual

Definition at line 116 of file CalWireSBND_module.cc.

  {
  }

Member Function Documentation

void caldata::CalWireSBND::beginJob

(

)

Definition at line 140 of file CalWireSBND_module.cc.

  {  
  }

void caldata::CalWireSBND::endJob

(

)

Definition at line 145 of file CalWireSBND_module.cc.

  {  
  }

void caldata::CalWireSBND::produce

(

art::Event &

evt

)

filter::ChannelFilter *chanFilt = new filter::ChannelFilter();

Definition at line 150 of file CalWireSBND_module.cc.

  {      
    // get the geometry
    art::ServiceHandle<geo::Geometry> geom;

    // get the FFT service to have access to the FFT size
    art::ServiceHandle<util::LArFFT> fFFT;
    int transformSize = fFFT->FFTSize();

    // Get signal shaping service.
    art::ServiceHandle<util::SignalShapingServiceSBND> sss;
    double DeconNorm = sss->GetDeconNorm();

    // make a collection of Wires
    std::unique_ptr<std::vector<recob::Wire> > wirecol(new std::vector<recob::Wire>);
    
    // ... and an association set     --Hec
    std::unique_ptr<art::Assns<raw::RawDigit,recob::Wire> > WireDigitAssn
      (new art::Assns<raw::RawDigit,recob::Wire>);
    
    // Read in the digit List object(s). 
    art::Handle< std::vector<raw::RawDigit> > digitVecHandle;
    if(fSpillName.size()>0) evt.getByLabel(fDigitModuleLabel, fSpillName, digitVecHandle);
    else evt.getByLabel(fDigitModuleLabel, digitVecHandle);

    if (!digitVecHandle->size())  return;
    mf::LogInfo("CalWireSBND") << "CalWireSBND:: digitVecHandle size is " << digitVecHandle->size();

    // Use the handle to get a particular (0th) element of collection.
    art::Ptr<raw::RawDigit> digitVec0(digitVecHandle, 0);
        
    unsigned int dataSize = digitVec0->Samples(); //size of raw data vectors


    if( (unsigned int)transformSize < dataSize){
      mf::LogWarning("CalWireSBND")<<"FFT size (" << transformSize << ") "
                                    << "is smaller than the data size (" << dataSize << ") "
                                    << "\nResizing the FFT now...";
      fFFT->ReinitializeFFT(dataSize,fFFT->FFTOptions(),fFFT->FFTFitBins());
      transformSize = fFFT->FFTSize();
      mf::LogWarning("CalWireSBND")<<"FFT size is now (" << transformSize << ") "
                                    << "and should be larger than the data size (" << dataSize << ")";
    }

    mf::LogInfo("CalWireSBND") << "Data size is " << dataSize << " and transform size is " << transformSize;

    if(fBaseSampleBins > 0 && dataSize % fBaseSampleBins != 0) {
      mf::LogError("CalWireSBND")<<"Set BaseSampleBins modulo dataSize= "<<dataSize;
    }

    uint32_t     channel(0); // channel number
    unsigned int bin(0);     // time bin loop variable
    
///    filter::ChannelFilter *chanFilt = new filter::ChannelFilter();  

    std::vector<float> holder;                // holds signal data
    std::vector<short> rawadc(transformSize);  // vector holding uncompressed adc values
    std::vector<TComplex> freqHolder(transformSize+1); // temporary frequency data
    
    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);

    // loop over all wires    
    wirecol->reserve(digitVecHandle->size());
    for(size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter){ // ++ move
      holder.clear();
      
      // get the reference to the current raw::RawDigit
      art::Ptr<raw::RawDigit> digitVec(digitVecHandle, rdIter);
      channel = digitVec->Channel();

      // skip bad channels
      //  if(!chanFilt->BadChannel(channel)) {
      if(true) {

        // resize and pad with zeros
        holder.resize(transformSize, 0.);
        
        // uncompress the data
        raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->Compression());
        
        // loop over all adc values and subtract the pedestal
        float pdstl = digitVec->GetPedestal();
        
        for(bin = 0; bin < dataSize; ++bin) 
          holder[bin]=(rawadc[bin]-pdstl);

        //fill the remaining bin with data
        for(bin = dataSize; bin < holder.size(); bin++){
          //  philosophy change - don't repeat data but instead fill extra space with zeros.
          //    not sure that one is better than the other.
          //      holder[bin] = (rawadc[bin-dataSize]-pdstl);
          holder[bin] = 0.0;
        }

        // Do deconvolution.
        sss->Deconvolute(clockData, channel, holder);
          for(bin = 0; bin < holder.size(); ++bin) holder[bin]=holder[bin]/DeconNorm;
      } // end if not a bad channel 
      
      holder.resize(dataSize,1e-5);

      // restore DC component through baseline subtraction
      if( fDoBaselineSub ) SubtractBaseline(holder);
      // more advanced, interpolation-based subtraction alg 
      // that uses the BaseSampleBins and BaseVarCut params
      if( fDoAdvBaselineSub ) SubtractBaselineAdv(holder);

      // Make a single ROI that spans the entire data size
      //RegionsOfInterest_t sparse_holder;
      //sparse_holder.add_range(0,holder.begin(),holder.end());
      CandidateROIVec candROIVec;
      fROITool->FindROIs( holder, channel, candROIVec);//calculates ROI and returns it to roiVec.
      //auto view = geom->View(channel);
      recob::Wire::RegionsOfInterest_t roiVec;

      //looping over roiVec to make a RegionOfInterest_t object.
      for(auto const& CandidateROI: candROIVec){
        size_t roiStart = CandidateROI.first;
        size_t roiStop = CandidateROI.second;
        std::vector<float> roiHolder;
        for(size_t i_holder=roiStart; i_holder<=roiStop; i_holder++){
          roiHolder.push_back(holder[i_holder]);
        }
        roiVec.add_range(roiStart, std::move(roiHolder));
      }
      wirecol->push_back(recob::WireCreator(std::move(roiVec),*digitVec).move());


      // add an association between the last object in wirecol--Hec
      // (that we just inserted) and digitVec
      if (!util::CreateAssn(*this, evt, *wirecol, digitVec, *WireDigitAssn, fSpillName)) {
        throw cet::exception("CalWireSBND")
          << "Can't associate wire #" << (wirecol->size() - 1)
          << " with raw digit #" << digitVec.key() << "\n";
      } // if failed to add association
    }


    if(wirecol->size() == 0)
      mf::LogWarning("CalWireSBND") << "No wires made for this event.";

    //--Hec if(fSpillName.size()>0)
    //--Hec   evt.put(std::move(wirecol), fSpillName);
    //--Hec else evt.put(std::move(wirecol));

    evt.put(std::move(wirecol), fSpillName);        //--Hec
    evt.put(std::move(WireDigitAssn), fSpillName);  //--Hec
    
   // delete chanFilt;
    return;
  }

void caldata::CalWireSBND::reconfigure

(

fhicl::ParameterSet const &

p

)

Definition at line 121 of file CalWireSBND_module.cc.

  {
    fDigitModuleLabel = p.get< std::string >("DigitModuleLabel", "daq");
    fDoBaselineSub    = p.get< bool >       ("DoBaselineSub");
    fDoAdvBaselineSub = p.get< bool >       ("DoAdvBaselineSub");
    fBaseSampleBins   = p.get< int >        ("BaseSampleBins");
    fBaseVarCut       = p.get< int >        ("BaseVarCut");
    
    fSpillName="";
    
    size_t pos = fDigitModuleLabel.find(":");
    if( pos!=std::string::npos ) {
      fSpillName = fDigitModuleLabel.substr( pos+1 );
      fDigitModuleLabel = fDigitModuleLabel.substr( 0, pos );
    }
    
  }

void caldata::CalWireSBND::SubtractBaseline ( std::vector< float > &  holder )

private

Definition at line 303 of file CalWireSBND_module.cc.

  {
    // Robust baseline calculation that effectively ignores outlier 
    // samples from large pulses:
    //   (1) fill a histogram with every sample's value,
    //   (2) find mode (bin with most entries),
    //   (3) calculate the mean along the entire waveform using
    //       only samples with values close to this mode.
    unsigned int bin(0);  
    float min = 0, max = 0;
    for(bin = 0; bin < holder.size(); bin++){
      if (holder[bin] > max) max = holder[bin];
      if (holder[bin] < min) min = holder[bin];
    }
    int nbin = max - min;
    if (nbin > 0) {
      TH1F h("h","h",nbin,min,max);
      for(bin = 0; bin < holder.size(); bin++) h.Fill(holder[bin]);
      float x_max = h.GetXaxis()->GetBinCenter(h.GetMaximumBin());
      float ped   = x_max;
      float sum   = 0;
      int ncount  = 0;
      for(bin = 0; bin < holder.size(); bin++){
        if( fabs(holder[bin]-x_max) < 2. ) {
          sum += holder[bin];
          ncount++; 
        }
      }
      if (ncount) ped = sum/ncount;
      for(bin = 0; bin < holder.size(); bin++) holder[bin] -= ped;
    }
  }

void caldata::CalWireSBND::SubtractBaselineAdv ( std::vector< float > &  holder )

private

Definition at line 336 of file CalWireSBND_module.cc.

  {
      // Subtract baseline using linear interpolation between regions defined
      // by the datasize and fBaseSampleBins

      // number of points to characterize the baseline
      unsigned short nBasePts = holder.size() / fBaseSampleBins;

      // the baseline offset vector
      std::vector<float> base;
      for(unsigned short ii = 0; ii < nBasePts; ++ii) base.push_back(0.);
      // find the average value in each region, using values that are
      // similar
      float fbins = fBaseSampleBins;
      unsigned short nfilld = 0;
      for(unsigned short ii = 0; ii < nBasePts; ++ii) {
        unsigned short loBin = ii * fBaseSampleBins;
        unsigned short hiBin = loBin + fBaseSampleBins;
        float ave = 0.;
        float sum = 0.;
        for(unsigned short bin = loBin; bin < hiBin; ++bin) {
          ave += holder[bin];
          sum += holder[bin] * holder[bin];
        } // jj
        ave = ave / fbins;
        float var = (sum - fbins * ave * ave) / (fbins - 1.);
        // Set the baseline for this region if the variance is small
        if(var < fBaseVarCut) {
          base[ii] = ave;
          ++nfilld;
        }
      } // ii
      // fill in any missing points if there aren't too many missing
      if(nfilld < nBasePts && nfilld > nBasePts / 2) {
        bool baseOK = true;
        // check the first region
        if(base[0] == 0) {
          unsigned short ii1 = 0;
          for(unsigned short ii = 1; ii < nBasePts; ++ii) {
            if(base[ii] != 0) {
              ii1 = ii;
              break;
            }
          } // ii
          unsigned short ii2 = 0;
          for(unsigned short ii = ii1 + 1; ii < nBasePts; ++ii) {
            if(base[ii] != 0) {
              ii2 = ii;
              break;
            }
          } // ii
          // failure
          if(ii2 > 0) {
            float slp = (base[ii2] - base[ii1]) / (float)(ii2 - ii1);
            base[0] = base[ii1] - slp * ii1;
          } else {
            baseOK = false;
          }
        } // base[0] == 0
        // check the last region
        if(baseOK && base[nBasePts] == 0) {
          unsigned short ii2 = 0;
          for(unsigned short ii = nBasePts - 1; ii > 0; --ii) {
            if(base[ii] != 0) {
              ii2 = ii;
              break;
            }
          } // ii
          baseOK = false; // assume the worst, hope for better
          unsigned short ii1 = 0;
          if (ii2 >= 1) {
            for(unsigned short ii = ii2 - 1; ii > 0; --ii) {
              if(base[ii] != 0) {
                ii1 = ii;
                baseOK = true;
                break;
              } // if base[ii]
            } // ii
          } // if ii2
          if (baseOK) {
            float slp = (base[ii2] - base[ii1]) / (float)(ii2 - ii1);
            base[nBasePts] = base[ii2] + slp * (nBasePts - ii2);
          }
        } // baseOK && base[nBasePts] == 0
        // now fill in any intermediate points
        for(unsigned short ii = 1; ii < nBasePts - 1; ++ii) {
          if(base[ii] == 0) {
            // find the next non-zero region
            for(unsigned short jj = ii + 1; jj < nBasePts; ++jj) {
              if(base[jj] != 0) {
                float slp = (base[jj] - base[ii - 1]) / (jj - ii + 1);
                base[ii] = base[ii - 1] + slp;
                break;
              }
            } // jj
          } // base[ii] == 0
        } // ii
      } // nfilld < nBasePts

      // interpolate and subtract
      float slp = (base[1] - base[0]) / (float)fBaseSampleBins;
      // bin offset to the origin (the center of the region)
      unsigned short bof = fBaseSampleBins / 2;
      unsigned short lastRegion = 0;
      for(unsigned short bin = 0; bin < holder.size(); ++bin) {
        // in a new region?
        unsigned short region = bin / fBaseSampleBins;
        if(region > lastRegion) {
          // update the slope and offset
          slp = (base[region] - base[lastRegion]) / (float)fBaseSampleBins;
          bof += fBaseSampleBins;
          lastRegion = region;
        }
        holder[bin] -= base[region] + (bin - bof) * slp;
      }
    }

Member Data Documentation

int caldata::CalWireSBND::fBaseSampleBins

private

bin grouping size in "interpolate" method

Definition at line 81 of file CalWireSBND_module.cc.

float caldata::CalWireSBND::fBaseVarCut

private

baseline variance cut used in "interpolate" method

Definition at line 82 of file CalWireSBND_module.cc.

std::string caldata::CalWireSBND::fDigitModuleLabel

private

module that made digits

Definition at line 84 of file CalWireSBND_module.cc.

bool caldata::CalWireSBND::fDoAdvBaselineSub

private

use interpolation-based baseline subtraction

Definition at line 80 of file CalWireSBND_module.cc.

bool caldata::CalWireSBND::fDoBaselineSub

private

subtract baseline to restore DC component post-deconvolution

Definition at line 79 of file CalWireSBND_module.cc.

std::unique_ptr<sbnd_tool::IROIFinder> caldata::CalWireSBND::fROITool

Definition at line 68 of file CalWireSBND_module.cc.

std::string caldata::CalWireSBND::fSpillName

private

nominal spill is an empty string it is set by the DigitModuleLabel ex.: "daq:preSpill" for prespill data

Definition at line 86 of file CalWireSBND_module.cc.

---

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

• CalWireSBND_module.cc