Inheritance diagram for caldata::RecoWireROIICARUS:

Public Member Functions
	RecoWireROIICARUS (fhicl::ParameterSet const &pset)

virtual	~RecoWireROIICARUS ()

void	produce (art::Event &evt)

void	beginJob ()

void	endJob ()

void	reconfigure (fhicl::ParameterSet const &p)

Private Member Functions
float	fixTheFreakingWaveform (const std::vector< float > &, raw::ChannelID_t, std::vector< float > &)

float	getTruncatedRMS (const std::vector< float > &) const

Private Attributes
std::string	fDigitModuleLabel
	module that made digits More...

std::string	fSpillName

unsigned short	fNoiseSource
	Used to determine ROI threshold. More...

int	fSaveWireWF
	Save recob::wire object waveforms. More...

size_t	fEventCount
	count of event processed More...

int	fMinAllowedChanStatus
	Don't consider channels with lower status. More...

float	fTruncRMSThreshold
	Calculate RMS up to this threshold... More...

float	fTruncRMSMinFraction
	or at least this fraction of time bins More...

bool	fOutputHistograms
	Output histograms? More...

std::vector< std::unique_ptr < icarus_tool::IROIFinder > >	fROIFinderVec
	ROI finders per plane. More...

std::unique_ptr < icarus_tool::IDeconvolution >	fDeconvolution

icarus_signal_processing::WaveformTools < float >	fWaveformTool

const geo::GeometryCore *	fGeometry = lar::providerFrom<geo::Geometry>()

std::vector< TH1F * >	fPedestalOffsetVec

std::vector< TH1F * >	fFullRMSVec

std::vector< TH1F * >	fTruncRMSVec

std::vector< TH1F * >	fNumTruncBinsVec

std::vector< TProfile * >	fPedByChanVec

std::vector< TProfile * >	fTruncRMSByChanVec

std::vector< TH1F * >	fNumROIsHistVec

std::vector< TH1F * >	fROILenHistVec

Detailed Description

Definition at line 61 of file RecoWireROIICARUS_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 115 of file RecoWireROIICARUS_module.cc.

                                                                   : EDProducer{pset}
 {
   this->reconfigure(pset);
 
   produces< std::vector<recob::Wire> >(fSpillName);
   produces<art::Assns<raw::RawDigit, recob::Wire>>(fSpillName);
 }

caldata::RecoWireROIICARUS::~RecoWireROIICARUS ( )

virtual

Definition at line 124 of file RecoWireROIICARUS_module.cc.

125 {

126 }

Member Function Documentation

void caldata::RecoWireROIICARUS::beginJob ( )

Definition at line 197 of file RecoWireROIICARUS_module.cc.

 {
     fEventCount = 0;
 } // beginJob

void caldata::RecoWireROIICARUS::endJob ( )

Definition at line 203 of file RecoWireROIICARUS_module.cc.

204 {

205 }

float caldata::RecoWireROIICARUS::fixTheFreakingWaveform	(	const std::vector< float > &	waveform,
		raw::ChannelID_t	channel,
		std::vector< float > &	fixedWaveform
	)

private

Definition at line 374 of file RecoWireROIICARUS_module.cc.

 {
     // Get the truncated mean and rms
     float fullRMS;
     float truncRMS;
     float truncMean;
     float nSig(2.0);  // make tight constraint
     int   nTrunc;
     int   range;
 
     fixedWaveform.resize(waveform.size());
     
     fWaveformTool.getPedestalCorrectedWaveform(waveform, fixedWaveform, nSig, truncMean, fullRMS, truncRMS, nTrunc, range);
     
     // Fill histograms
     if (fOutputHistograms)
     {
         std::vector<geo::WireID> wids = fGeometry->ChannelToWire(channel);
     
         // Recover plane and wire in the plane
         size_t plane = wids[0].Plane;
         size_t wire  = wids[0].Wire;
     
         fPedestalOffsetVec[plane]->Fill(truncMean,1.);
         fFullRMSVec[plane]->Fill(fullRMS, 1.);
         fTruncRMSVec[plane]->Fill(truncRMS, 1.);
         fNumTruncBinsVec[plane]->Fill(nTrunc, 1.);
         fPedByChanVec[plane]->Fill(wire, truncMean, 1.);
         fTruncRMSByChanVec[plane]->Fill(wire, truncRMS, 1.);
     }
     
     return truncRMS;
 }

float caldata::RecoWireROIICARUS::getTruncatedRMS ( const std::vector< float > & waveform ) const

private

Definition at line 352 of file RecoWireROIICARUS_module.cc.

 {
     // do rms calculation - the old fashioned way and over all adc values
     std::vector<float> locWaveform = waveform;
     
     // sort in ascending order so we can truncate the sume
     std::sort(locWaveform.begin(), locWaveform.end(),[](const auto& left, const auto& right){return std::fabs(left) < std::fabs(right);});
     
     float threshold = fTruncRMSThreshold;
     
     std::vector<float>::iterator threshItr = std::find_if(locWaveform.begin(),locWaveform.end(),[threshold](const auto& val){return std::fabs(val) > threshold;});
     
     int minNumBins = std::max(int(fTruncRMSMinFraction * locWaveform.size()),int(std::distance(locWaveform.begin(),threshItr)));
     
     // Get the truncated sum
     float truncRms = std::inner_product(locWaveform.begin(), locWaveform.begin() + minNumBins, locWaveform.begin(), 0.);
     
     truncRms = std::sqrt(std::max(0.,truncRms / double(minNumBins)));
     
     return truncRms;
 }

void caldata::RecoWireROIICARUS::produce ( art::Event & evt )

Definition at line 208 of file RecoWireROIICARUS_module.cc.

 {
     //get pedestal conditions
     const lariov::DetPedestalProvider& pedestalRetrievalAlg = art::ServiceHandle<lariov::DetPedestalService>()->GetPedestalProvider();
     
     // make a collection of Wires
     std::unique_ptr<std::vector<recob::Wire> > wirecol(new std::vector<recob::Wire>);
     // ... and an association set
     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())
     {
         evt.put(std::move(wirecol), fSpillName);
         evt.put(std::move(WireDigitAssn), fSpillName);
         return;
     }
     
     raw::ChannelID_t channel = raw::InvalidChannelID; // channel number
     
     const lariov::ChannelStatusProvider& chanFilt = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
     
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
     double const samplingRate = sampling_rate(clockData);
     // loop over all wires
     wirecol->reserve(digitVecHandle->size());
     for(size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter)
     {
         // vector that will be moved into the Wire object
         recob::Wire::RegionsOfInterest_t ROIVec;
 
         // get the reference to the current raw::RawDigit
         art::Ptr<raw::RawDigit> digitVec(digitVecHandle, rdIter);
         channel = digitVec->Channel();
       
         // The following test is meant to be temporary until the "correct" solution is implemented
         if (!chanFilt.IsPresent(channel)) continue;
 
         // Testing an idea about rejecting channels
         if (digitVec->GetPedestal() < 0.) continue;
 
         float pedestal = 0.;
         
         // skip bad channels
         if( chanFilt.Status(channel) >= fMinAllowedChanStatus)
         {
             size_t dataSize = digitVec->Samples();
             
             // Recover the plane info
             std::vector<geo::WireID> wids    = fGeometry->ChannelToWire(channel);
             const geo::PlaneID&      planeID = wids[0].planeID();
 
             // vector holding uncompressed adc values
             std::vector<short> rawadc(dataSize);
             
             // uncompress the data
             raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->Compression());
             
             // loop over all adc values and subtract the pedestal
             // When we have a pedestal database, can provide the digit timestamp as the third argument of GetPedestalMean
             pedestal = pedestalRetrievalAlg.PedMean(channel);
             
             // Get the pedestal subtracted data, centered in the deconvolution vector
             std::vector<float> rawAdcLessPedVec(dataSize);
             
             std::transform(rawadc.begin(),rawadc.end(),rawAdcLessPedVec.begin(),std::bind(std::minus<short>(),std::placeholders::_1,pedestal));
             
             // It seems there are deviations from the pedestal when using wirecell for noise filtering
             float raw_noise = fixTheFreakingWaveform(rawAdcLessPedVec, channel, rawAdcLessPedVec);
             
             // Recover a measure of the noise on the channel for use in the ROI finder
             //float raw_noise = getTruncatedRMS(rawAdcLessPedVec);
             
             // Try smoothing the input waveform
 //            std::vector<float> rawAdcSmoothVec;
 //            fWaveformTool->medianSmooth(rawAdcLessPedVec,rawAdcSmoothVec);
             
             // vector of candidate ROI begin and end bins
             icarus_tool::IROIFinder::CandidateROIVec candRoiVec;
 
             // Now find the candidate ROI's
             fROIFinderVec.at(planeID.Plane)->FindROIs(rawAdcLessPedVec, channel, fEventCount, raw_noise, candRoiVec);
             
             // Do the deconvolution
             fDeconvolution->Deconvolve(rawAdcLessPedVec, samplingRate, channel, candRoiVec, ROIVec);
             
             // Make some histograms?
             if (fOutputHistograms)
             {
                 // First up, determine what kind of wire we have
                 std::vector<geo::WireID> wids    = fGeometry->ChannelToWire(channel);
                 const geo::PlaneID&      planeID = wids[0].planeID();
                 
                 fNumROIsHistVec.at(planeID.Plane)->Fill(candRoiVec.size(), 1.);
                 
                 for(const auto& pair : candRoiVec)
                     fROILenHistVec.at(planeID.Plane)->Fill(pair.second-pair.first, 1.);
             }
         } // end if not a bad channel
 
         // create the new wire directly in wirecol
         wirecol->push_back(recob::WireCreator(std::move(ROIVec),*digitVec).move());
 
         // add an association between the last object in wirecol
         // (that we just inserted) and digitVec
         if (!util::CreateAssn(*this, evt, *wirecol, digitVec, *WireDigitAssn, fSpillName))
         {
             throw art::Exception(art::errors::ProductRegistrationFailure)
                 << "Can't associate wire #" << (wirecol->size() - 1)
                 << " with raw digit #" << digitVec.key();
         } // if failed to add association
         //  DumpWire(wirecol->back()); // for debugging
     }
 
     if(wirecol->size() == 0)
       mf::LogWarning("RecoWireROIICARUS") << "No wires made for this event.";
 
     //Make Histogram of recob::wire objects from Signal() vector
     // get access to the TFile service
     if ( fSaveWireWF ){
         art::ServiceHandle<art::TFileService> tfs;
         for (size_t wireN = 0; wireN < wirecol->size(); wireN++){
             std::vector<float> sigTMP = wirecol->at(wireN).Signal();
             TH1D* fWire = tfs->make<TH1D>(Form("Noise_Evt%04zu_N%04zu",fEventCount,wireN), ";Noise (ADC);",
                                       sigTMP.size(),-0.5,sigTMP.size()-0.5);
             for (size_t tick = 0; tick < sigTMP.size(); tick++){
                 fWire->SetBinContent(tick+1, sigTMP.at(tick) );
             }
         }
     }
     
     evt.put(std::move(wirecol), fSpillName);
     evt.put(std::move(WireDigitAssn), fSpillName);
 
     fEventCount++;
 
     return;
 } // produce

void caldata::RecoWireROIICARUS::reconfigure ( fhicl::ParameterSet const & p )

Definition at line 129 of file RecoWireROIICARUS_module.cc.

 {
     // Recover the vector of fhicl parameters for the ROI tools
     const fhicl::ParameterSet& roiFinderTools = pset.get<fhicl::ParameterSet>("ROIFinderToolVec");
     
     fROIFinderVec.resize(roiFinderTools.get_pset_names().size());
     
     for(const std::string& roiFinderTool : roiFinderTools.get_pset_names())
     {
         const fhicl::ParameterSet& roiFinderToolParamSet = roiFinderTools.get<fhicl::ParameterSet>(roiFinderTool);
         size_t                     planeIdx              = roiFinderToolParamSet.get<size_t>("Plane");
         
         fROIFinderVec.at(planeIdx) = art::make_tool<icarus_tool::IROIFinder>(roiFinderToolParamSet);
     }
     
     std::sort(fROIFinderVec.begin(),fROIFinderVec.end(),[](const auto& left,const auto& right){return left->plane() < right->plane();});
 
     fDeconvolution = art::make_tool<icarus_tool::IDeconvolution>(pset.get<fhicl::ParameterSet>("Deconvolution"));
 
     fDigitModuleLabel           = pset.get< std::string >   ("DigitModuleLabel", "daq");
     fNoiseSource                = pset.get< unsigned short >("NoiseSource",          3);
     fSaveWireWF                 = pset.get< int >           ("SaveWireWF"             );
     fMinAllowedChanStatus       = pset.get< int >           ("MinAllowedChannelStatus");
     fTruncRMSThreshold          = pset.get< float >         ("TruncRMSThreshold",    6.);
     fTruncRMSMinFraction        = pset.get< float >         ("TruncRMSMinFraction", 0.6);
     fOutputHistograms           = pset.get< bool  >         ("OutputHistograms",   true);
     
     fSpillName.clear();
     
     size_t pos = fDigitModuleLabel.find(":");
     if( pos!=std::string::npos )
     {
         fSpillName = fDigitModuleLabel.substr( pos+1 );
         fDigitModuleLabel = fDigitModuleLabel.substr( 0, pos );
     }
     
     if (fOutputHistograms)
     {
         // Access ART's TFileService, which will handle creating and writing
         // histograms and n-tuples for us.
         art::ServiceHandle<art::TFileService> tfs;
     
         fPedestalOffsetVec.resize(3);
         fFullRMSVec.resize(3);
         fTruncRMSVec.resize(3);
         fNumTruncBinsVec.resize(3);
         fPedByChanVec.resize(3);
         fTruncRMSByChanVec.resize(3);
         fNumROIsHistVec.resize(3);
         fROILenHistVec.resize(3);
     
         for(size_t planeIdx = 0; planeIdx < 3; planeIdx++)
         {
             fPedestalOffsetVec[planeIdx] = tfs->make<TH1F>(    Form("PedPlane_%02zu",planeIdx),            ";Pedestal Offset (ADC);", 100, -5., 5.);
             fFullRMSVec[planeIdx]        = tfs->make<TH1F>(    Form("RMSFPlane_%02zu",planeIdx),           "Full RMS;RMS (ADC);", 100, 0., 10.);
             fTruncRMSVec[planeIdx]       = tfs->make<TH1F>(    Form("RMSTPlane_%02zu",planeIdx),           "Truncated RMS;RMS (ADC);", 100, 0., 10.);
             fNumTruncBinsVec[planeIdx]   = tfs->make<TH1F>(    Form("NTruncBins_%02zu",planeIdx),          ";# bins",     640, 0., 6400.);
             fPedByChanVec[planeIdx]      = tfs->make<TProfile>(Form("PedByWirePlane_%02zu",planeIdx),      ";Wire#", fGeometry->Nwires(planeIdx), 0., fGeometry->Nwires(planeIdx), -5., 5.);
             fTruncRMSByChanVec[planeIdx] = tfs->make<TProfile>(Form("TruncRMSByWirePlane_%02zu",planeIdx), ";Wire#", fGeometry->Nwires(planeIdx), 0., fGeometry->Nwires(planeIdx),  0., 10.);
             fNumROIsHistVec[planeIdx]    = tfs->make<TH1F>(    Form("NROISplane_%02zu",planeIdx),          ";# ROIs;",   100, 0, 100);
             fROILenHistVec[planeIdx]     = tfs->make<TH1F>(    Form("ROISIZEplane_%02zu",planeIdx),        ";ROI size;", 500, 0, 500);
         }
     }
     
     return;
 }