Inheritance diagram for hit::FFTHitFinder:

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

Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDProducer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept

	Producer (fhicl::ParameterSet const &)

	Producer (Producer const &)=delete

	Producer (Producer &&)=delete

Producer &	operator= (Producer const &)=delete

Producer &	operator= (Producer &&)=delete

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

Public Member Functions inherited from art::Modifier
	~Modifier () noexcept

	Modifier ()

	Modifier (Modifier const &)=delete

	Modifier (Modifier &&)=delete

Modifier &	operator= (Modifier const &)=delete

Modifier &	operator= (Modifier &&)=delete

Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept

	ModuleBase ()

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

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

Private Attributes
std::string	fCalDataModuleLabel

double	fMinSigInd
	Induction signal height threshold. More...

double	fMinSigCol
	Collection signal height threshold. More...

double	fIndWidth
	Initial width for induction fit. More...

double	fColWidth
	Initial width for collection fit. More...

double	fIndMinWidth
	Minimum induction hit width. More...

double	fColMinWidth
	Minimum collection hit width. More...

int	fMaxMultiHit
	maximum hits for multi fit More...

int	fAreaMethod
	Type of area calculation. More...

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

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer

using	WorkerType = WorkerT< EDProducer >

Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >

Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >

Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)

Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 38 of file FFTHitFinder_module.cc.

Constructor & Destructor Documentation

hit::FFTHitFinder::FFTHitFinder ( fhicl::ParameterSet const & pset )

explicit

Definition at line 61 of file FFTHitFinder_module.cc.

     : EDProducer{pset}
   {
     fCalDataModuleLabel = pset.get< std::string  >("CalDataModuleLabel");
     fMinSigInd          = pset.get< double       >("MinSigInd");
     fMinSigCol          = pset.get< double       >("MinSigCol");
     fIndWidth           = pset.get< double       >("IndWidth");
     fColWidth           = pset.get< double       >("ColWidth");
     fIndMinWidth        = pset.get< double       >("IndMinWidth");
     fColMinWidth        = pset.get< double       >("ColMinWidth");
     fMaxMultiHit        = pset.get< int          >("MaxMultiHit");
     fAreaMethod         = pset.get< int          >("AreaMethod");
     fAreaNorms          = pset.get< std::vector< double > >("AreaNorms");
 
     // let HitCollectionCreator declare that we are going to produce
     // hits and associations with wires and raw digits
     // (with no particular product label)
     recob::HitCollectionCreator::declare_products(producesCollector());
   }

Member Function Documentation

void hit::FFTHitFinder::produce ( art::Event & evt )

overrideprivatevirtual

Todo:

just get the integral from the fit for totSig

Todo:: need to have a disambiguation algorithm somewhere in here

Todo:

multiplicity and local_index have to be determined

Implements art::EDProducer.

Definition at line 86 of file FFTHitFinder_module.cc.

   {
 
     // this object contains the hit collection
     // and its associations to wires and raw digits:
     recob::HitCollectionCreator hcol(evt);
 
     // Read in the wire List object(s).
     art::Handle< std::vector<recob::Wire> > wireVecHandle;
     evt.getByLabel(fCalDataModuleLabel,wireVecHandle);
     art::ServiceHandle<geo::Geometry const> geom;
 
     // also get the raw digits associated with wires
     art::FindOneP<raw::RawDigit> WireToRawDigits
       (wireVecHandle, evt, fCalDataModuleLabel);
 
     std::vector<int> startTimes;             // stores time of 1st local minimum
     std::vector<int> maxTimes;               // stores time of local maximum
     std::vector<int> endTimes;               // stores time of 2nd local minimum
     int time               = 0;              // current time bin
     int minTimeHolder      = 0;              // current start time
     raw::ChannelID_t channel = raw::InvalidChannelID; // channel number
     bool maxFound          = false;          // Flag for whether a peak > threshold has been found
     double threshold       = 0.;             // minimum signal size for id'ing a hit
     double fitWidth        = 0.;             // hit fit width initial value
     double minWidth        = 0.;             // minimum hit width
     std::string eqn        = "gaus(0)";      // string for equation to fit
     geo::SigType_t sigType = geo::kInduction;// type of plane we are looking at
     std::stringstream numConv;
 
 
     //loop over wires
     for(unsigned int wireIter = 0; wireIter < wireVecHandle->size(); wireIter++) {
 
       art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
       startTimes.clear();
       maxTimes.clear();
       endTimes.clear();
       std::vector<float> signal(wire->Signal());
       std::vector<float>::iterator timeIter;   // iterator for time bins
       time          = 0;
       minTimeHolder = 0;
       maxFound      = false;
       channel       = wire->Channel();
       sigType       = geom->SignalType(channel);
 
       //Set the appropriate signal widths and thresholds
       if(sigType == geo::kInduction){
         threshold     = fMinSigInd;
         fitWidth      = fIndWidth;
         minWidth      = fIndMinWidth;
       }
       else if(sigType == geo::kCollection){
         threshold = fMinSigCol;
         fitWidth  = fColWidth;
         minWidth  = fColMinWidth;
       }
       // loop over signal
       for(timeIter = signal.begin(); timeIter+2 < signal.end(); timeIter++){
         //test if timeIter+1 is a local minimum
         if(*timeIter > *(timeIter+1) && *(timeIter+1) < *(timeIter+2)){
           //only add points if already found a local max above threshold.
           if(maxFound) {
             endTimes.push_back(time+1);
             maxFound = false;
             //keep these in case new hit starts right away
             minTimeHolder = time+2;
           }
           else minTimeHolder = time+1;
         }
         //if not a minimum, test if we are at a local maximum
         //if so, and the max value is above threshold, add it and proceed.
         else if(*timeIter < *(timeIter+1) &&
                 *(timeIter+1) > *(timeIter+2) &&
                 *(timeIter+1) > threshold){
           maxFound = true;
           maxTimes.push_back(time+1);
           startTimes.push_back(minTimeHolder);
         }
         time++;
       }//end loop over signal vec
 
 
       //if no inflection found before end, but peak found add end point
       while(maxTimes.size()>endTimes.size())
         endTimes.push_back(signal.size()-1);
       if(startTimes.size() == 0) continue;
 
       //All code below does the fitting, adding of hits
       //to the hit vector and when all wires are complete
       //saving them
       double totSig(0); // stores the total hit signal
       double startT(0); // stores the start time
       double endT(0);   // stores the end time
       int numHits(0);   // number of consecutive hits being fitted
       int size(0);      // size of data vector for fit
       int hitIndex(0);  // index of current hit in sequence
       double amplitude(0), position(0), width(0);  //fit parameters
       double amplitudeErr(0), positionErr(0), widthErr(0);  //fit errors
       double goodnessOfFit(0), chargeErr(0);  //Chi2/NDF and error on charge
       double minPeakHeight(0);  //lowest peak height in multi-hit fit
 
       //stores gaussian paramters first index is the hit number
       //the second refers to height, position, and width respectively
       std::vector<double>  hitSig;
 
       //add found hits to hit vector
       while(hitIndex < (signed)startTimes.size()) {
 
         startT = endT = 0;
         numHits = 1;
         minPeakHeight = signal[maxTimes[hitIndex]];
 
         //consider adding pulse to group of consecutive hits if:
         //1 less than max consecutive hits
         //2 we are not at the last point in the signal vector
         //3 the height of the dip between the two is greater than threshold/2
         //4 and there is no gap between them
         while(numHits < fMaxMultiHit &&
               numHits+hitIndex < (signed)endTimes.size() &&
               signal[endTimes[hitIndex+numHits-1]] >threshold/2.0 &&
               startTimes[hitIndex+numHits] - endTimes[hitIndex+numHits-1] < 2){
 
           if(signal[maxTimes[hitIndex+numHits]] < minPeakHeight)
             minPeakHeight = signal[maxTimes[hitIndex+numHits]];
 
           ++numHits;
         }
 
         //finds the first point > 1/2 the smallest peak
         startT = startTimes[hitIndex];
 
         while(signal[(int)startT] < minPeakHeight/2.0) ++startT;
 
         //finds the first point from the end > 1/2 the smallest peak
         endT = endTimes[hitIndex+numHits-1];
 
         while(signal[(int)endT] <minPeakHeight/2.0) --endT;
         size = (int)(endT-startT);
         TH1D hitSignal("hitSignal","",size,startT,endT);
         for(int i = (int)startT; i < (int)endT; ++i)
           hitSignal.Fill(i,signal[i]);
 
         //build the TFormula
         eqn = "gaus(0)";
 
         for(int i = 3; i < numHits*3; i+=3) {
           eqn.append("+gaus(");
           numConv.str("");
           numConv << i;
           eqn.append(numConv.str());
           eqn.append(")");
         }
 
         TF1 gSum("gSum",eqn.c_str(),0,size);
 
         if(numHits > 1) {
           TArrayD data(numHits*numHits);
           TVectorD amps(numHits);
           for(int i = 0; i < numHits; ++i) {
             amps[i] = signal[maxTimes[hitIndex+i]];
             for(int j = 0; j < numHits;j++)
               data[i+numHits*j] = TMath::Gaus(maxTimes[hitIndex+j],
                                               maxTimes[hitIndex+i],
                                               fitWidth);
           }//end loop over hits
 
           //This section uses a linear approximation in order to get an
           //initial value of the individual hit amplitudes
           try{
             TMatrixD h(numHits,numHits);
             h.Use(numHits,numHits,data.GetArray());
             TDecompSVD a(h);
             a.Solve(amps);
           }
           catch(...){
             mf::LogInfo("FFTHitFinder")<<"TDcompSVD failed";
             hitIndex += numHits;
             continue;
           }
 
           for(int i = 0; i < numHits; ++i) {
             //if the approximation makes a peak vanish
             //set initial height as average of threshold and
             //raw peak height
             if(amps[i] > 0 ) amplitude = amps[i];
             else amplitude = 0.5*(threshold+signal[maxTimes[hitIndex+i]]);
             gSum.SetParameter(3*i,amplitude);
             gSum.SetParameter(1+3*i, maxTimes[hitIndex+i]);
             gSum.SetParameter(2+3*i, fitWidth);
             gSum.SetParLimits(3*i, 0.0, 3.0*amplitude);
             gSum.SetParLimits(1+3*i, startT , endT);
             gSum.SetParLimits(2+3*i, 0.0, 10.0*fitWidth);
           }//end loop over hits
         }//end if numHits > 1
         else {
           gSum.SetParameters(signal[maxTimes[hitIndex]],maxTimes[hitIndex],fitWidth);
           gSum.SetParLimits(0,0.0,1.5*signal[maxTimes[hitIndex]]);
           gSum.SetParLimits(1, startT , endT);
           gSum.SetParLimits(2,0.0,10.0*fitWidth);
         }
 
         /// \todo - just get the integral from the fit for totSig
         hitSignal.Fit(&gSum,"QNRW","", startT, endT);
         for(int hitNumber = 0; hitNumber < numHits; ++hitNumber) {
           totSig = 0;
           if(gSum.GetParameter(3*hitNumber)   > threshold/2.0 &&
              gSum.GetParameter(3*hitNumber+2) > minWidth) {
             amplitude     = gSum.GetParameter(3*hitNumber);
             position      = gSum.GetParameter(3*hitNumber+1);
             width         = gSum.GetParameter(3*hitNumber+2);
             amplitudeErr  = gSum.GetParError(3*hitNumber);
             positionErr   = gSum.GetParError(3*hitNumber+1);
             widthErr      = gSum.GetParError(3*hitNumber+2);
             goodnessOfFit = gSum.GetChisquare()/(double)gSum.GetNDF();
        int DoF = gSum.GetNDF();
 
             //estimate error from area of Gaussian
             chargeErr = std::sqrt(TMath::Pi())*(amplitudeErr*width+widthErr*amplitude);
 
             hitSig.resize(size);
 
             for(int sigPos = 0; sigPos < size; ++sigPos){
               hitSig[sigPos] = amplitude*TMath::Gaus(sigPos+startT,position, width);
               totSig += hitSig[(int)sigPos];
             }
 
             if(fAreaMethod)
               totSig = std::sqrt(2*TMath::Pi())*amplitude*width/fAreaNorms[(size_t)sigType];
 
             // get the WireID for this hit
             std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
             ///\todo need to have a disambiguation algorithm somewhere in here
             // for now, just take the first option returned from ChannelToWire
             geo::WireID wid = wids[0];
 
             // make the hit
             recob::HitCreator hit(
               *wire,          // wire
               wid,            // wireID
               (int) startT,   // start_tick
               (int) endT,     // end_tick
               width,          // rms
               position,       // peak_time
               positionErr,    // sigma_peak_time
               amplitude,      // peak_amplitude
               amplitudeErr,   // sigma_peak_amplitude
               totSig,         // hit_integral
               chargeErr,      // hit_sigma_integral
               std::accumulate // summedADC
                 (signal.begin() + (int) startT, signal.begin() + (int) endT, 0.),
               1,              // multiplicity
               -1,             // local_index
                               /// \todo - multiplicity and local_index have to be determined
               goodnessOfFit,  // goodness_of_fit
               DoF             // dof
               );
 
             // get the object associated with the original hit
             art::Ptr<raw::RawDigit> rawdigits = WireToRawDigits.at(wireIter);
 
             hcol.emplace_back(hit.move(), wire, rawdigits);
           }//end if over threshold
         }//end loop over hits
         hitIndex += numHits;
       } // end while on hitIndex<(signed)startTimes.size()
 
     } // while on Wires
 
     // put the hit collection and associations into the event
     hcol.put_into(evt);
 
   } // End of produce()

Member Data Documentation

int hit::FFTHitFinder::fAreaMethod

private

Type of area calculation.

Definition at line 55 of file FFTHitFinder_module.cc.

std::vector<double> hit::FFTHitFinder::fAreaNorms

private

factors for converting area to same units as peak height

Definition at line 56 of file FFTHitFinder_module.cc.

std::string hit::FFTHitFinder::fCalDataModuleLabel

private

Definition at line 47 of file FFTHitFinder_module.cc.

double hit::FFTHitFinder::fColMinWidth

private

Minimum collection hit width.

Definition at line 53 of file FFTHitFinder_module.cc.

double hit::FFTHitFinder::fColWidth

private

Initial width for collection fit.

Definition at line 51 of file FFTHitFinder_module.cc.

double hit::FFTHitFinder::fIndMinWidth

private

Minimum induction hit width.

Definition at line 52 of file FFTHitFinder_module.cc.

double hit::FFTHitFinder::fIndWidth

private

Initial width for induction fit.

Definition at line 50 of file FFTHitFinder_module.cc.

int hit::FFTHitFinder::fMaxMultiHit

private

maximum hits for multi fit

Definition at line 54 of file FFTHitFinder_module.cc.

double hit::FFTHitFinder::fMinSigCol

private

Collection signal height threshold.

Definition at line 49 of file FFTHitFinder_module.cc.

double hit::FFTHitFinder::fMinSigInd

private

Induction signal height threshold.

Definition at line 48 of file FFTHitFinder_module.cc.

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

larreco/larreco/HitFinder/FFTHitFinder_module.cc

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation