Inheritance diagram for hit::GausHitFinder:

Public Member Functions
	GausHitFinder (fhicl::ParameterSet const &pset, art::ProcessingFrame const &)

Private Member Functions
void	produce (art::Event &evt, art::ProcessingFrame const &) override

void	beginJob (art::ProcessingFrame const &) override

std::vector< double >	FillOutHitParameterVector (const std::vector< double > &input)

Private Attributes
const bool	fFilterHits

const bool	fFillHists

const std::string	fCalDataModuleLabel

const std::string	fAllHitsInstanceName

const std::vector< int >	fLongMaxHitsVec
	Maximum number hits on a really long pulse train. More...

const std::vector< int >	fLongPulseWidthVec
	Sets width of hits used to describe long pulses. More...

const size_t	fMaxMultiHit
	maximum hits for multi fit More...

const int	fAreaMethod
	Type of area calculation. More...

const std::vector< double >	fAreaNormsVec
	factors for converting area to same units as peak height More...

const double	fChi2NDF

const std::vector< float >	fPulseHeightCuts
	maximum Chisquared / NDF allowed for a hit to be saved More...

const std::vector< float >	fPulseWidthCuts

const std::vector< float >	fPulseRatioCuts

std::atomic< size_t >	fEventCount {0}

std::vector< std::unique_ptr < reco_tool::ICandidateHitFinder > >	fHitFinderToolVec
	For finding candidate hits. More...

std::unique_ptr < reco_tool::IPeakFitter >	fPeakFitterTool
	Perform fit to candidate peaks. More...

std::unique_ptr< HitFilterAlg >	fHitFilterAlg
	algorithm used to filter out noise hits More...

TH1F *	fFirstChi2

TH1F *	fChi2

Detailed Description

Definition at line 62 of file GausHitFinder_module.cc.

Constructor & Destructor Documentation

hit::GausHitFinder::GausHitFinder	(	fhicl::ParameterSet const &	pset,
		art::ProcessingFrame const &
	)

explicit

Definition at line 109 of file GausHitFinder_module.cc.

     : SharedProducer{pset}
     , fFilterHits(pset.get<bool>("FilterHits", false))
     , fFillHists(pset.get<bool>("FillHists", false))
     , fCalDataModuleLabel(pset.get<std::string>("CalDataModuleLabel"))
     , fAllHitsInstanceName(pset.get<std::string>("AllHitsInstanceName", ""))
     , fLongMaxHitsVec(pset.get<std::vector<int>>("LongMaxHits", std::vector<int>() = {25, 25, 25}))

Member Function Documentation

void hit::GausHitFinder::beginJob ( art::ProcessingFrame const & )

overrideprivate

Definition at line 201 of file GausHitFinder_module.cc.

   {
     // get access to the TFile service
     art::ServiceHandle<art::TFileService const> tfs;
 
     // ======================================
     // === Hit Information for Histograms ===
     if (fFillHists) {
       fFirstChi2 = tfs->make<TH1F>("fFirstChi2", "#chi^{2}", 10000, 0, 5000);
       fChi2 = tfs->make<TH1F>("fChi2", "#chi^{2}", 10000, 0, 5000);
     }
   }

std::vector< double > hit::GausHitFinder::FillOutHitParameterVector ( const std::vector< double > & input )

private

Definition at line 178 of file GausHitFinder_module.cc.

   {
     if (input.size() == 0)
       throw std::runtime_error(
         "GausHitFinder::FillOutHitParameterVector ERROR! Input config vector has zero size.");
 
     std::vector<double> output;
     art::ServiceHandle<geo::Geometry const> geom;
     const unsigned int N_PLANES = geom->Nplanes();
 
     if (input.size() == 1)
       output.resize(N_PLANES, input[0]);
     else if (input.size() == N_PLANES)
       output = input;
     else
       throw std::runtime_error("GausHitFinder::FillOutHitParameterVector ERROR! Input config "
                                "vector size !=1 and !=N_PLANES.");
     return output;
   }

void hit::GausHitFinder::produce	(	art::Event &	evt,
		art::ProcessingFrame const &
	)

overrideprivate

Definition at line 219 of file GausHitFinder_module.cc.

   {
     unsigned int count = fEventCount.fetch_add(1);
     //==================================================================================================
 
     TH1::AddDirectory(kFALSE);
 
     // Instantiate and Reset a stop watch
     //TStopwatch StopWatch;
     //StopWatch.Reset();
 
     // ################################
     // ### Calling Geometry service ###
     // ################################
     art::ServiceHandle<geo::Geometry const> geom;
 
     // ###############################################
     // ### Making a ptr vector to put on the event ###
     // ###############################################
     // this contains the hit collection
     // and its associations to wires and raw digits
     recob::HitCollectionCreator allHitCol(evt, fAllHitsInstanceName, true, false);
 
     // Handle the filtered hits collection...
     recob::HitCollectionCreator hcol(evt, "", true, false);
     recob::HitCollectionCreator* filteredHitCol = 0;
 
     if (fFilterHits) filteredHitCol = &hcol;
 
     //store in a thread safe way
     struct hitstruct {
       recob::Hit hit_tbb;
       art::Ptr<recob::Wire> wire_tbb;
     };
 
     tbb::concurrent_vector<hitstruct> hitstruct_vec;
     tbb::concurrent_vector<hitstruct> filthitstruct_vec;
 
     //    if (fAllHitsInstanceName != "") filteredHitCol = &hcol;
 
     // ##########################################
     // ### Reading in the Wire List object(s) ###
     // ##########################################
     art::Handle<std::vector<recob::Wire>> wireVecHandle;
     evt.getByLabel(fCalDataModuleLabel, wireVecHandle);
 
     //#################################################
     //###    Set the charge determination method    ###
     //### Default is to compute the normalized area ###
     //#################################################
     std::function<double(double, double, double, double, int, int)> chargeFunc =
       [](double peakMean, double peakAmp, double peakWidth, double areaNorm, int low, int hi) {
         return std::sqrt(2 * TMath::Pi()) * peakAmp * peakWidth / areaNorm;
       };
 
     //##############################################
     //### Alternative is to integrate over pulse ###
     //##############################################
     if (fAreaMethod == 0)
       chargeFunc =
         [](double peakMean, double peakAmp, double peakWidth, double areaNorm, int low, int hi) {
           double charge(0);
           for (int sigPos = low; sigPos < hi; sigPos++)
             charge += peakAmp * TMath::Gaus(sigPos, peakMean, peakWidth);
           return charge;
         };
 
     //##############################
     //### Looping over the wires ###
     //##############################
     //for(size_t wireIter = 0; wireIter < wireVecHandle->size(); wireIter++)
     //{
     tbb::parallel_for(
       static_cast<std::size_t>(0),
       wireVecHandle->size(),
       [&](size_t& wireIter) {
         // ####################################
         // ### Getting this particular wire ###
         // ####################################
         art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
 
         // --- Setting Channel Number and Signal type ---
 
         raw::ChannelID_t channel = wire->Channel();
 
         // get the WireID for this hit
         std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
         // for now, just take the first option returned from ChannelToWire
         geo::WireID wid = wids[0];
         // We need to know the plane to look up parameters
         geo::PlaneID::PlaneID_t plane = wid.Plane;
 
         // ----------------------------------------------------------
         // -- Setting the appropriate signal widths and thresholds --
         // --    for the right plane.      --
         // ----------------------------------------------------------
 
         // #################################################
         // ### Set up to loop over ROI's for this wire   ###
         // #################################################
         const recob::Wire::RegionsOfInterest_t& signalROI = wire->SignalROI();
 
         // for (const auto& range : signalROI.get_ranges()) {
         tbb::parallel_for(static_cast<std::size_t>(0),signalROI.n_ranges(),
                           [&](size_t& rangeIter){
           const auto& range = signalROI.range(rangeIter);
           // ROI start time
           raw::TDCtick_t roiFirstBinTick = range.begin_index();
 
           // ###########################################################
           // ### Scan the waveform and find candidate peaks + merge  ###
           // ###########################################################
 
           reco_tool::ICandidateHitFinder::HitCandidateVec hitCandidateVec;
           reco_tool::ICandidateHitFinder::MergeHitCandidateVec mergedCandidateHitVec;
 
           fHitFinderToolVec.at(plane)->findHitCandidates(range, 0, channel, count, hitCandidateVec);
           fHitFinderToolVec.at(plane)->MergeHitCandidates(
             range, hitCandidateVec, mergedCandidateHitVec);
 
           // #######################################################
           // ### Lets loop over the pulses we found on this wire ###
           // #######################################################
 
           for (auto& mergedCands : mergedCandidateHitVec) {
             int startT = mergedCands.front().startTick;
             int endT = mergedCands.back().stopTick;
 
             // ### Putting in a protection in case things went wrong ###
             // ### In the end, this primarily catches the case where ###
             // ### a fake pulse is at the start of the ROI           ###
             if (endT - startT < 5) continue;
 
             // #######################################################
             // ### Clearing the parameter vector for the new pulse ###
             // #######################################################
 
             // === Setting the number of Gaussians to try ===
             int nGausForFit = mergedCands.size();
 
             // ##################################################
             // ### Calling the function for fitting Gaussians ###
             // ##################################################
             double chi2PerNDF(0.);
             int NDF(1);
             /*stand alone
                 reco_tool::IPeakFitter::PeakParamsVec peakParamsVec(nGausForFit);
                 */
             reco_tool::IPeakFitter::PeakParamsVec peakParamsVec;
 
             // #######################################################
             // ### If # requested Gaussians is too large then punt ###
             // #######################################################
             if (mergedCands.size() <= fMaxMultiHit) {
               fPeakFitterTool->findPeakParameters(
                 range.data(), mergedCands, peakParamsVec, chi2PerNDF, NDF);
 
               // If the chi2 is infinite then there is a real problem so we bail
               if (!(chi2PerNDF < std::numeric_limits<double>::infinity())) {
                 chi2PerNDF = 2. * fChi2NDF;
                 NDF = 2;
               }
 
               if (fFillHists) fFirstChi2->Fill(chi2PerNDF);
             }
 
             // #######################################################
             // ### If too large then force alternate solution      ###
             // ### - Make n hits from pulse train where n will     ###
             // ###   depend on the fhicl parameter fLongPulseWidth ###
             // ### Also do this if chi^2 is too large              ###
             // #######################################################
             if (mergedCands.size() > fMaxMultiHit || nGausForFit * chi2PerNDF > fChi2NDF) {
               int longPulseWidth = fLongPulseWidthVec.at(plane);
               int nHitsThisPulse = (endT - startT) / longPulseWidth;
 
               if (nHitsThisPulse > fLongMaxHitsVec.at(plane)) {
                 nHitsThisPulse = fLongMaxHitsVec.at(plane);
                 longPulseWidth = (endT - startT) / nHitsThisPulse;
               }
 
               if (nHitsThisPulse * longPulseWidth < endT - startT) nHitsThisPulse++;
 
               int firstTick = startT;
               int lastTick = std::min(firstTick + longPulseWidth, endT);
 
               peakParamsVec.clear();
               nGausForFit = nHitsThisPulse;
               NDF = 1.;
               chi2PerNDF = chi2PerNDF > fChi2NDF ? chi2PerNDF : -1.;
 
               for (int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) {
                 // This hit parameters
                 double sumADC =
                   std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.);
                 double peakSigma = (lastTick - firstTick) / 3.; // Set the width...
                 double peakAmp = 0.3989 * sumADC / peakSigma;   // Use gaussian formulation
                 double peakMean = (firstTick + lastTick) / 2.;
 
                 // Store hit params
                 reco_tool::IPeakFitter::PeakFitParams_t peakParams;
 
                 peakParams.peakCenter = peakMean;
                 peakParams.peakCenterError = 0.1 * peakMean;
                 peakParams.peakSigma = peakSigma;
                 peakParams.peakSigmaError = 0.1 * peakSigma;
                 peakParams.peakAmplitude = peakAmp;
                 peakParams.peakAmplitudeError = 0.1 * peakAmp;
 
                 peakParamsVec.push_back(peakParams);
 
                 // set for next loop
                 firstTick = lastTick;
                 lastTick = std::min(lastTick + longPulseWidth, endT);
               }
             }
 
             // #######################################################
             // ### Loop through returned peaks and make recob hits ###
             // #######################################################
 
             int numHits(0);
 
             // Make a container for what will be the filtered collection
             std::vector<recob::Hit> filteredHitVec;
 
             for (const auto& peakParams : peakParamsVec) {
               // Extract values for this hit
               float peakAmp = peakParams.peakAmplitude;
               float peakMean = peakParams.peakCenter;
               float peakWidth = peakParams.peakSigma;
 
               // Place one bit of protection here
               if (std::isnan(peakAmp)) {
                 std::cout << "**** hit peak amplitude is a nan! Channel: " << channel
                           << ", start tick: " << startT << std::endl;
                 continue;
               }
 
               // Extract errors
               float peakAmpErr = peakParams.peakAmplitudeError;
               float peakMeanErr = peakParams.peakCenterError;
               float peakWidthErr = peakParams.peakSigmaError;
 
               // ### Charge ###
               float charge =
                 chargeFunc(peakMean, peakAmp, peakWidth, fAreaNormsVec[plane], startT, endT);
               ;
               float chargeErr =
                 std::sqrt(TMath::Pi()) * (peakAmpErr * peakWidthErr + peakWidthErr * peakAmpErr);
 
               // ### limits for getting sums
               std::vector<float>::const_iterator sumStartItr = range.begin() + startT;
               std::vector<float>::const_iterator sumEndItr = range.begin() + endT;
 
               // ### Sum of ADC counts
               double sumADC = std::accumulate(sumStartItr, sumEndItr, 0.);
 
               // ok, now create the hit
               recob::HitCreator hitcreator(
                 *wire,                      // wire reference
                 wid,                        // wire ID
                 startT + roiFirstBinTick,   // start_tick TODO check
                 endT + roiFirstBinTick,     // end_tick TODO check
                 peakWidth,                  // rms
                 peakMean + roiFirstBinTick, // peak_time
                 peakMeanErr,                // sigma_peak_time
                 peakAmp,                    // peak_amplitude
                 peakAmpErr,                 // sigma_peak_amplitude
                 charge,                     // hit_integral
                 chargeErr,                  // hit_sigma_integral
                 sumADC,                     // summedADC FIXME
                 nGausForFit,                // multiplicity
                 numHits,                    // local_index TODO check that the order is correct
                 chi2PerNDF,                 // goodness_of_fit
                 NDF                         // dof
               );
 
               if (filteredHitCol) filteredHitVec.push_back(hitcreator.copy());
 
               const recob::Hit hit(hitcreator.move());
 
               // This loop will store ALL hits
               hitstruct tmp{std::move(hit), wire};
               hitstruct_vec.push_back(std::move(tmp));
 
               numHits++;
             } // <---End loop over gaussians
 
             // Should we filter hits?
             if (filteredHitCol && !filteredHitVec.empty()) {
               // #######################################################################
               // Is all this sorting really necessary?  Would it be faster to just loop
               // through the hits and perform simple cuts on amplitude and width on a
               // hit-by-hit basis, either here in the module (using fPulseHeightCuts and
               // fPulseWidthCuts) or in HitFilterAlg?
               // #######################################################################
 
               // Sort in ascending peak height
               std::sort(filteredHitVec.begin(),
                         filteredHitVec.end(),
                         [](const auto& left, const auto& right) {
                           return left.PeakAmplitude() > right.PeakAmplitude();
                         });
 
               // Reject if the first hit fails the PH/wid cuts
               if (filteredHitVec.front().PeakAmplitude() < fPulseHeightCuts.at(plane) ||
                   filteredHitVec.front().RMS() < fPulseWidthCuts.at(plane))
                 filteredHitVec.clear();
 
               // Now check other hits in the snippet
               if (filteredHitVec.size() > 1) {
                 // The largest pulse height will now be at the front...
                 float largestPH = filteredHitVec.front().PeakAmplitude();
 
                 // Find where the pulse heights drop below threshold
                 float threshold(fPulseRatioCuts.at(plane));
 
                 std::vector<recob::Hit>::iterator smallHitItr = std::find_if(
                   filteredHitVec.begin(),
                   filteredHitVec.end(),
                   [largestPH, threshold](const auto& hit) {
                     return hit.PeakAmplitude() < 8. && hit.PeakAmplitude() / largestPH < threshold;
                   });
 
                 // Shrink to fit
                 if (smallHitItr != filteredHitVec.end())
                   filteredHitVec.resize(std::distance(filteredHitVec.begin(), smallHitItr));
 
                 // Resort in time order
                 std::sort(filteredHitVec.begin(),
                           filteredHitVec.end(),
                           [](const auto& left, const auto& right) {
                             return left.PeakTime() < right.PeakTime();
                           });
               }
 
               // Copy the hits we want to keep to the filtered hit collection
               for (const auto& filteredHit : filteredHitVec)
                 if (!fHitFilterAlg || fHitFilterAlg->IsGoodHit(filteredHit)) {
                   hitstruct tmp{std::move(filteredHit), wire};
                   filthitstruct_vec.push_back(std::move(tmp));
                 }
 
               if (fFillHists) fChi2->Fill(chi2PerNDF);
             }
           } //<---End loop over merged candidate hits
 
         } //<---End looping over ROI's
         );//end tbb parallel for
 
       }   //<---End looping over all the wires
     );    //end tbb parallel for
 
     for (size_t i = 0; i < hitstruct_vec.size(); i++) {
       allHitCol.emplace_back(hitstruct_vec[i].hit_tbb, hitstruct_vec[i].wire_tbb);
     }
 
     for (size_t j = 0; j < filthitstruct_vec.size(); j++) {
       filteredHitCol->emplace_back(filthitstruct_vec[j].hit_tbb, filthitstruct_vec[j].wire_tbb);
     }
 
     //==================================================================================================
     // End of the event -- move the hit collection and the associations into the event
 
     if (filteredHitCol) {
 
       // If we filtered hits but no instance name was
       // specified for the "all hits" collection, then
       // only save the filtered hits to the event
       if (fAllHitsInstanceName == "") {
         filteredHitCol->put_into(evt);
 
         // otherwise, save both
       }
       else {
         filteredHitCol->put_into(evt);
         allHitCol.put_into(evt);
       }
     }
     else {
       allHitCol.put_into(evt);
     }
 
     // Keep track of events processed
     //fEventCount++;
 
   } // End of produce()

Member Data Documentation

const std::string hit::GausHitFinder::fAllHitsInstanceName

private

Definition at line 76 of file GausHitFinder_module.cc.

const int hit::GausHitFinder::fAreaMethod

private

Type of area calculation.

Definition at line 82 of file GausHitFinder_module.cc.

const std::vector<double> hit::GausHitFinder::fAreaNormsVec

private

factors for converting area to same units as peak height

Definition at line 84 of file GausHitFinder_module.cc.

const std::string hit::GausHitFinder::fCalDataModuleLabel

private

Definition at line 75 of file GausHitFinder_module.cc.

TH1F* hit::GausHitFinder::fChi2

private

Definition at line 103 of file GausHitFinder_module.cc.

const double hit::GausHitFinder::fChi2NDF

private

Definition at line 85 of file GausHitFinder_module.cc.

std::atomic<size_t> hit::GausHitFinder::fEventCount {0}

private

Definition at line 91 of file GausHitFinder_module.cc.

const bool hit::GausHitFinder::fFillHists

private

Definition at line 73 of file GausHitFinder_module.cc.

const bool hit::GausHitFinder::fFilterHits

private

Definition at line 72 of file GausHitFinder_module.cc.

TH1F* hit::GausHitFinder::fFirstChi2

private

Definition at line 102 of file GausHitFinder_module.cc.

std::unique_ptr<HitFilterAlg> hit::GausHitFinder::fHitFilterAlg

private

algorithm used to filter out noise hits

Definition at line 99 of file GausHitFinder_module.cc.

std::vector<std::unique_ptr<reco_tool::ICandidateHitFinder> > hit::GausHitFinder::fHitFinderToolVec

private

For finding candidate hits.

Definition at line 95 of file GausHitFinder_module.cc.

const std::vector<int> hit::GausHitFinder::fLongMaxHitsVec

private

Maximum number hits on a really long pulse train.

Definition at line 78 of file GausHitFinder_module.cc.

const std::vector<int> hit::GausHitFinder::fLongPulseWidthVec

private

Sets width of hits used to describe long pulses.

Definition at line 79 of file GausHitFinder_module.cc.

const size_t hit::GausHitFinder::fMaxMultiHit

private

maximum hits for multi fit

Definition at line 81 of file GausHitFinder_module.cc.

std::unique_ptr<reco_tool::IPeakFitter> hit::GausHitFinder::fPeakFitterTool

private

Perform fit to candidate peaks.

Definition at line 97 of file GausHitFinder_module.cc.

const std::vector<float> hit::GausHitFinder::fPulseHeightCuts

private

maximum Chisquared / NDF allowed for a hit to be saved

Definition at line 87 of file GausHitFinder_module.cc.

const std::vector<float> hit::GausHitFinder::fPulseRatioCuts

private

Definition at line 89 of file GausHitFinder_module.cc.

const std::vector<float> hit::GausHitFinder::fPulseWidthCuts

private

Definition at line 88 of file GausHitFinder_module.cc.

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

GausHitFinder_module.cc

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation