Public Member Functions | Private Member Functions | Private Attributes | List of all members
hit::FFTHitFinder Class Reference
Inheritance diagram for hit::FFTHitFinder:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

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
 
Produceroperator= (Producer const &)=delete
 
Produceroperator= (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
 
Modifieroperator= (Modifier const &)=delete
 
Modifieroperator= (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
ConsumesCollectorconsumesCollector ()
 
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.

62  : EDProducer{pset}
63  {
64  fCalDataModuleLabel = pset.get< std::string >("CalDataModuleLabel");
65  fMinSigInd = pset.get< double >("MinSigInd");
66  fMinSigCol = pset.get< double >("MinSigCol");
67  fIndWidth = pset.get< double >("IndWidth");
68  fColWidth = pset.get< double >("ColWidth");
69  fIndMinWidth = pset.get< double >("IndMinWidth");
70  fColMinWidth = pset.get< double >("ColMinWidth");
71  fMaxMultiHit = pset.get< int >("MaxMultiHit");
72  fAreaMethod = pset.get< int >("AreaMethod");
73  fAreaNorms = pset.get< std::vector< double > >("AreaNorms");
74 
75  // let HitCollectionCreator declare that we are going to produce
76  // hits and associations with wires and raw digits
77  // (with no particular product label)
79  }
double fIndMinWidth
Minimum induction hit width.
std::string string
Definition: nybbler.cc:12
double fMinSigInd
Induction signal height threshold.
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.h:20
double fIndWidth
Initial width for induction fit.
int fMaxMultiHit
maximum hits for multi fit
double fColWidth
Initial width for collection fit.
static void declare_products(art::ProducesCollector &collector, std::string instance_name="", bool doWireAssns=true, bool doRawDigitAssns=true)
Declares the hit products we are going to fill.
Definition: HitCreator.cxx:248
std::vector< double > fAreaNorms
factors for converting area to same units as peak height
std::string fCalDataModuleLabel
double fColMinWidth
Minimum collection hit width.
ProducesCollector & producesCollector() noexcept
double fMinSigCol
Collection signal height threshold.
int fAreaMethod
Type of area calculation.

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.

87  {
88 
89  // this object contains the hit collection
90  // and its associations to wires and raw digits:
92 
93  // Read in the wire List object(s).
95  evt.getByLabel(fCalDataModuleLabel,wireVecHandle);
97 
98  // also get the raw digits associated with wires
99  art::FindOneP<raw::RawDigit> WireToRawDigits
100  (wireVecHandle, evt, fCalDataModuleLabel);
101 
102  std::vector<int> startTimes; // stores time of 1st local minimum
103  std::vector<int> maxTimes; // stores time of local maximum
104  std::vector<int> endTimes; // stores time of 2nd local minimum
105  int time = 0; // current time bin
106  int minTimeHolder = 0; // current start time
107  raw::ChannelID_t channel = raw::InvalidChannelID; // channel number
108  bool maxFound = false; // Flag for whether a peak > threshold has been found
109  double threshold = 0.; // minimum signal size for id'ing a hit
110  double fitWidth = 0.; // hit fit width initial value
111  double minWidth = 0.; // minimum hit width
112  std::string eqn = "gaus(0)"; // string for equation to fit
113  geo::SigType_t sigType = geo::kInduction;// type of plane we are looking at
114  std::stringstream numConv;
115 
116 
117  //loop over wires
118  for(unsigned int wireIter = 0; wireIter < wireVecHandle->size(); wireIter++) {
119 
120  art::Ptr<recob::Wire> wire(wireVecHandle, wireIter);
121  startTimes.clear();
122  maxTimes.clear();
123  endTimes.clear();
124  std::vector<float> signal(wire->Signal());
125  std::vector<float>::iterator timeIter; // iterator for time bins
126  time = 0;
127  minTimeHolder = 0;
128  maxFound = false;
129  channel = wire->Channel();
130  sigType = geom->SignalType(channel);
131 
132  //Set the appropriate signal widths and thresholds
133  if(sigType == geo::kInduction){
134  threshold = fMinSigInd;
135  fitWidth = fIndWidth;
136  minWidth = fIndMinWidth;
137  }
138  else if(sigType == geo::kCollection){
139  threshold = fMinSigCol;
140  fitWidth = fColWidth;
141  minWidth = fColMinWidth;
142  }
143  // loop over signal
144  for(timeIter = signal.begin(); timeIter+2 < signal.end(); timeIter++){
145  //test if timeIter+1 is a local minimum
146  if(*timeIter > *(timeIter+1) && *(timeIter+1) < *(timeIter+2)){
147  //only add points if already found a local max above threshold.
148  if(maxFound) {
149  endTimes.push_back(time+1);
150  maxFound = false;
151  //keep these in case new hit starts right away
152  minTimeHolder = time+2;
153  }
154  else minTimeHolder = time+1;
155  }
156  //if not a minimum, test if we are at a local maximum
157  //if so, and the max value is above threshold, add it and proceed.
158  else if(*timeIter < *(timeIter+1) &&
159  *(timeIter+1) > *(timeIter+2) &&
160  *(timeIter+1) > threshold){
161  maxFound = true;
162  maxTimes.push_back(time+1);
163  startTimes.push_back(minTimeHolder);
164  }
165  time++;
166  }//end loop over signal vec
167 
168 
169  //if no inflection found before end, but peak found add end point
170  while(maxTimes.size()>endTimes.size())
171  endTimes.push_back(signal.size()-1);
172  if(startTimes.size() == 0) continue;
173 
174  //All code below does the fitting, adding of hits
175  //to the hit vector and when all wires are complete
176  //saving them
177  double totSig(0); // stores the total hit signal
178  double startT(0); // stores the start time
179  double endT(0); // stores the end time
180  int numHits(0); // number of consecutive hits being fitted
181  int size(0); // size of data vector for fit
182  int hitIndex(0); // index of current hit in sequence
183  double amplitude(0), position(0), width(0); //fit parameters
184  double amplitudeErr(0), positionErr(0), widthErr(0); //fit errors
185  double goodnessOfFit(0), chargeErr(0); //Chi2/NDF and error on charge
186  double minPeakHeight(0); //lowest peak height in multi-hit fit
187 
188  //stores gaussian paramters first index is the hit number
189  //the second refers to height, position, and width respectively
190  std::vector<double> hitSig;
191 
192  //add found hits to hit vector
193  while(hitIndex < (signed)startTimes.size()) {
194 
195  startT = endT = 0;
196  numHits = 1;
197  minPeakHeight = signal[maxTimes[hitIndex]];
198 
199  //consider adding pulse to group of consecutive hits if:
200  //1 less than max consecutive hits
201  //2 we are not at the last point in the signal vector
202  //3 the height of the dip between the two is greater than threshold/2
203  //4 and there is no gap between them
204  while(numHits < fMaxMultiHit &&
205  numHits+hitIndex < (signed)endTimes.size() &&
206  signal[endTimes[hitIndex+numHits-1]] >threshold/2.0 &&
207  startTimes[hitIndex+numHits] - endTimes[hitIndex+numHits-1] < 2){
208 
209  if(signal[maxTimes[hitIndex+numHits]] < minPeakHeight)
210  minPeakHeight = signal[maxTimes[hitIndex+numHits]];
211 
212  ++numHits;
213  }
214 
215  //finds the first point > 1/2 the smallest peak
216  startT = startTimes[hitIndex];
217 
218  while(signal[(int)startT] < minPeakHeight/2.0) ++startT;
219 
220  //finds the first point from the end > 1/2 the smallest peak
221  endT = endTimes[hitIndex+numHits-1];
222 
223  while(signal[(int)endT] <minPeakHeight/2.0) --endT;
224  size = (int)(endT-startT);
225  TH1D hitSignal("hitSignal","",size,startT,endT);
226  for(int i = (int)startT; i < (int)endT; ++i)
227  hitSignal.Fill(i,signal[i]);
228 
229  //build the TFormula
230  eqn = "gaus(0)";
231 
232  for(int i = 3; i < numHits*3; i+=3) {
233  eqn.append("+gaus(");
234  numConv.str("");
235  numConv << i;
236  eqn.append(numConv.str());
237  eqn.append(")");
238  }
239 
240  TF1 gSum("gSum",eqn.c_str(),0,size);
241 
242  if(numHits > 1) {
243  TArrayD data(numHits*numHits);
244  TVectorD amps(numHits);
245  for(int i = 0; i < numHits; ++i) {
246  amps[i] = signal[maxTimes[hitIndex+i]];
247  for(int j = 0; j < numHits;j++)
248  data[i+numHits*j] = TMath::Gaus(maxTimes[hitIndex+j],
249  maxTimes[hitIndex+i],
250  fitWidth);
251  }//end loop over hits
252 
253  //This section uses a linear approximation in order to get an
254  //initial value of the individual hit amplitudes
255  try{
256  TMatrixD h(numHits,numHits);
257  h.Use(numHits,numHits,data.GetArray());
258  TDecompSVD a(h);
259  a.Solve(amps);
260  }
261  catch(...){
262  mf::LogInfo("FFTHitFinder")<<"TDcompSVD failed";
263  hitIndex += numHits;
264  continue;
265  }
266 
267  for(int i = 0; i < numHits; ++i) {
268  //if the approximation makes a peak vanish
269  //set initial height as average of threshold and
270  //raw peak height
271  if(amps[i] > 0 ) amplitude = amps[i];
272  else amplitude = 0.5*(threshold+signal[maxTimes[hitIndex+i]]);
273  gSum.SetParameter(3*i,amplitude);
274  gSum.SetParameter(1+3*i, maxTimes[hitIndex+i]);
275  gSum.SetParameter(2+3*i, fitWidth);
276  gSum.SetParLimits(3*i, 0.0, 3.0*amplitude);
277  gSum.SetParLimits(1+3*i, startT , endT);
278  gSum.SetParLimits(2+3*i, 0.0, 10.0*fitWidth);
279  }//end loop over hits
280  }//end if numHits > 1
281  else {
282  gSum.SetParameters(signal[maxTimes[hitIndex]],maxTimes[hitIndex],fitWidth);
283  gSum.SetParLimits(0,0.0,1.5*signal[maxTimes[hitIndex]]);
284  gSum.SetParLimits(1, startT , endT);
285  gSum.SetParLimits(2,0.0,10.0*fitWidth);
286  }
287 
288  /// \todo - just get the integral from the fit for totSig
289  hitSignal.Fit(&gSum,"QNRW","", startT, endT);
290  for(int hitNumber = 0; hitNumber < numHits; ++hitNumber) {
291  totSig = 0;
292  if(gSum.GetParameter(3*hitNumber) > threshold/2.0 &&
293  gSum.GetParameter(3*hitNumber+2) > minWidth) {
294  amplitude = gSum.GetParameter(3*hitNumber);
295  position = gSum.GetParameter(3*hitNumber+1);
296  width = gSum.GetParameter(3*hitNumber+2);
297  amplitudeErr = gSum.GetParError(3*hitNumber);
298  positionErr = gSum.GetParError(3*hitNumber+1);
299  widthErr = gSum.GetParError(3*hitNumber+2);
300  goodnessOfFit = gSum.GetChisquare()/(double)gSum.GetNDF();
301  int DoF = gSum.GetNDF();
302 
303  //estimate error from area of Gaussian
304  chargeErr = std::sqrt(TMath::Pi())*(amplitudeErr*width+widthErr*amplitude);
305 
306  hitSig.resize(size);
307 
308  for(int sigPos = 0; sigPos < size; ++sigPos){
309  hitSig[sigPos] = amplitude*TMath::Gaus(sigPos+startT,position, width);
310  totSig += hitSig[(int)sigPos];
311  }
312 
313  if(fAreaMethod)
314  totSig = std::sqrt(2*TMath::Pi())*amplitude*width/fAreaNorms[(size_t)sigType];
315 
316  // get the WireID for this hit
317  std::vector<geo::WireID> wids = geom->ChannelToWire(channel);
318  ///\todo need to have a disambiguation algorithm somewhere in here
319  // for now, just take the first option returned from ChannelToWire
320  geo::WireID wid = wids[0];
321 
322  // make the hit
324  *wire, // wire
325  wid, // wireID
326  (int) startT, // start_tick
327  (int) endT, // end_tick
328  width, // rms
329  position, // peak_time
330  positionErr, // sigma_peak_time
331  amplitude, // peak_amplitude
332  amplitudeErr, // sigma_peak_amplitude
333  totSig, // hit_integral
334  chargeErr, // hit_sigma_integral
335  std::accumulate // summedADC
336  (signal.begin() + (int) startT, signal.begin() + (int) endT, 0.),
337  1, // multiplicity
338  -1, // local_index
339  /// \todo - multiplicity and local_index have to be determined
340  goodnessOfFit, // goodness_of_fit
341  DoF // dof
342  );
343 
344  // get the object associated with the original hit
345  art::Ptr<raw::RawDigit> rawdigits = WireToRawDigits.at(wireIter);
346 
347  hcol.emplace_back(hit.move(), wire, rawdigits);
348  }//end if over threshold
349  }//end loop over hits
350  hitIndex += numHits;
351  } // end while on hitIndex<(signed)startTimes.size()
352 
353  } // while on Wires
354 
355  // put the hit collection and associations into the event
356  hcol.put_into(evt);
357 
358  } // End of produce()
double fIndMinWidth
Minimum induction hit width.
intermediate_table::iterator iterator
std::string string
Definition: nybbler.cc:12
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
double fMinSigInd
Induction signal height threshold.
double fIndWidth
Initial width for induction fit.
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
int fMaxMultiHit
maximum hits for multi fit
double fColWidth
Initial width for collection fit.
uint8_t channel
Definition: CRTFragment.hh:201
SigType_t SignalType(geo::PlaneID const &pid) const
Returns the type of signal on the channels of specified TPC plane.
std::vector< double > fAreaNorms
factors for converting area to same units as peak height
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:83
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:633
std::string fCalDataModuleLabel
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:32
A class handling a collection of hits and its associations.
Definition: HitCreator.h:508
Signal from induction planes.
Definition: geo_types.h:145
enum geo::_plane_sigtype SigType_t
const double a
double fColMinWidth
Minimum collection hit width.
Detector simulation of raw signals on wires.
double fMinSigCol
Collection signal height threshold.
int fAreaMethod
Type of area calculation.
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
Signal from collection planes.
Definition: geo_types.h:146

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: