LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
RawHitFinder_module.cc
Go to the documentation of this file.
1 //HIT FINDER THAT RUNS ON RAW SIGNALS INSTEAD OF DECONVOLUTED ONES.
3 //WRITTEN INITIALLY FOR DUNE 35T ONLINE FILTER.
4 //abooth@fnal.gov, mstancar@fnal.gov.
6 
7 // C/C++ standard libraries
8 #include <numeric> // std::accumulate
9 #include <string>
10 
11 //Framework
19 #include "fhiclcpp/ParameterSet.h"
21 
22 //LArSoft
25 #include "lardataobj/RawData/raw.h"
26 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
27 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
28 
29 //LArSoft From FFT
30 
33 
34 namespace hit {
35 
36  class RawHitFinder : public art::EDProducer {
37 
38  public:
39  explicit RawHitFinder(fhicl::ParameterSet const& pset);
40 
41  private:
42  void produce(art::Event& evt) override;
43 
44  unsigned int fDataSize; //SIZE OF RAW DATA ON ONE WIRE.
45  art::InputTag fDigitModuleLabel; //MODULE THAT MADE DIGITS.
46  std::string fSpillName; //NOMINAL SPILL IS AN EMPTY STRING.
47 
48  //FFT COPIED VARIABLES.
49  std::string fCalDataModuleLabel;
50  std::string fHitLabelName;
51  double fMinSigInd; //INDUCTION SIGNAL HEIGHT THRESHOLD.
52  double fMinSigCol; //COLLECTION SIGNAL HEIGHT THRESHOLD.
53  double fIndWidth; //INITIAL WIDTH FOR COLLECTION FIT.
54  double fColWidth; //INITIAL WIDTH FOR COLLECTION FIT.
55  double fIndMinWidth; //MINIMUM INDUCTION HIT WIDTH.
56  double fColMinWidth; //MINIMUM COLLECTION HIT WIDTH.
57  int fMaxMultiHit; //MAXIMUM HITS FOR MULTIFIT.
58  int fAreaMethod; //TYPE OF AREA CALCULATION.
59  std::vector<double> fAreaNorms; //FACTORS FOR CONVERTING AREA TO SAME UNIT AS PEAK HEIGHT.
60  bool fUncompressWithPed; //OPTION TO UNCOMPRESS WITH PEDESTAL.
61  bool fSkipInd; //OPTION TO SKIP INDUCTION PLANES
62  double fIncludeMoreTail; //FRACTION OF THE HIT WIDTH INCLUDED BELOW THRESHOLD IN
63  // THE CALCULATION OF CHARGE. COLLECTION PLANE ONLY
64  int fColMinWindow; // Minimum length of integration window for charge in ticks
65  int fIndCutoff; //MAX WIDTH FOR EARLY SIDE OF INDUCTION HIT IN TICKS
66 
67  }; // class RawHitFinder
68 
69  //-------------------------------------------------
71  {
72  fDigitModuleLabel = pset.get<art::InputTag>("DigitModuleLabel", "daq");
73  fCalDataModuleLabel = pset.get<std::string>("CalDataModuleLabel");
74  fMinSigInd = pset.get<double>("MinSigInd");
75  fMinSigCol = pset.get<double>("MinSigCol");
76  fIncludeMoreTail = pset.get<double>("IncludeMoreTail", 0.);
77  fIndWidth = pset.get<int>("IndWidth", 20);
78  fColWidth = pset.get<double>("ColWidth");
79  fIndMinWidth = pset.get<double>("IndMinWidth");
80  fColMinWidth = pset.get<double>("ColMinWidth", 0.);
81  fMaxMultiHit = pset.get<int>("MaxMultiHit");
82  fAreaMethod = pset.get<int>("AreaMethod");
83  fAreaNorms = pset.get<std::vector<double>>("AreaNorms");
84  fUncompressWithPed = pset.get<bool>("UncompressWithPed", true);
85  fSkipInd = pset.get<bool>("SkipInd", false);
86  fColMinWindow = pset.get<int>("ColMinWindow", 0);
87  fIndCutoff = pset.get<int>("IndCutoff", 20);
88  mf::LogInfo("RawHitFinder_module") << "fDigitModuleLabel: " << fDigitModuleLabel << std::endl;
89 
90  //LET HITCOLLECTIONCREATOR DECLARE THAT WE ARE GOING TO PRODUCE
91  //HITS AND ASSOCIATIONS TO RAW DIGITS BUT NOT ASSOCIATIONS TO WIRES
92  //(WITH NO PARTICULAR PRODUCT LABEL).
94  /*instance_name*/ "",
95  /*doWireAssns*/ false,
96  /*doRawDigitAssns*/ true);
97  }
98 
99  //-------------------------------------------------
101  {
102  //READ IN THE DIGIT LIST OBJECTS(S).
104 
105  bool retVal = evt.getByLabel(fDigitModuleLabel, digitVecHandle);
106  if (retVal == true)
107  mf::LogInfo("RawHitFinder_module")
108  << "I got fDigitModuleLabel: " << fDigitModuleLabel << std::endl;
109  else
110  mf::LogWarning("RawHitFinder_module")
111  << "Could not get fDigitModuleLabel: " << fDigitModuleLabel << std::endl;
112 
113  std::vector<float> holder; //HOLDS SIGNAL DATA.
114  std::vector<short> rawadc; //UNCOMPRESSED ADC VALUES.
115 
116  // ###############################################
117  // ### Making a ptr vector to put on the event ###
118  // ###############################################
119  // THIS CONTAINS THE HIT COLLECTION AND ITS ASSOCIATIONS TO WIRES AND RAW DIGITS.
120  recob::HitCollectionCreator hcol(evt, false /* doWireAssns */, true /* doRawDigitAssns */);
121 
122  std::vector<float> startTimes; //STORES TIME OF WINDOW START.
123  std::vector<float> maxTimes; //STORES TIME OF LOCAL MAXIMUM.
124  std::vector<float> endTimes; //STORES TIME OF WINDOW END.
125  std::vector<float> peakHeight; //STORES ADC COUNT AT THE MAXIMUM.
126  std::vector<float> hitrms; //STORES CHARGE WEIGHTED RMS OF TIME ACROSS THE HIT.
127  std::vector<double> charge; //STORES THE TOTAL CHARGE ASSOCIATED WITH THE HIT.
128 
129  uint32_t channel = 0; //CHANNEL NUMBER.
130  double threshold = 0; //MINIMUM SIGNAL SIZE FOR ID'ING A HIT.
131  double totSig = 0;
132  double myrms = 0;
133  double mynorm = 0;
135  std::stringstream numConv;
136 
137  hcol.reserve(digitVecHandle->size());
138  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
139  for (size_t rdIter = 0; rdIter < digitVecHandle->size(); ++rdIter) {
140  holder.clear();
141 
142  //GET THE REFERENCE TO THE CURRENT raw::RawDigit.
143  art::Ptr<raw::RawDigit> digitVec(digitVecHandle, rdIter);
144  channel = digitVec->Channel();
145  fDataSize = digitVec->Samples();
146 
147  rawadc.resize(fDataSize);
148  holder.resize(fDataSize);
149 
150  //UNCOMPRESS THE DATA.
151  if (fUncompressWithPed) {
152  int pedestal = (int)digitVec->GetPedestal();
153  raw::Uncompress(digitVec->ADCs(), rawadc, pedestal, digitVec->Compression());
154  }
155  else {
156  raw::Uncompress(digitVec->ADCs(), rawadc, digitVec->Compression());
157  }
158 
159  //GET THE LIST OF BAD CHANNELS.
160  lariov::ChannelStatusProvider const& channelStatus =
162 
163  lariov::ChannelStatusProvider::ChannelSet_t const BadChannels = channelStatus.BadChannels();
164 
165  for (unsigned int bin = 0; bin < fDataSize; ++bin) {
166  holder[bin] = (rawadc[bin] - digitVec->GetPedestal());
167  }
168 
169  sigType = wireReadoutGeom.SignalType(channel);
170 
171  peakHeight.clear();
172  endTimes.clear();
173  startTimes.clear();
174  maxTimes.clear();
175  charge.clear();
176  hitrms.clear();
177 
178  bool channelSwitch = false;
179 
180  for (auto it = BadChannels.begin(); it != BadChannels.end(); it++) {
181  if (channel == *it) {
182  channelSwitch = true;
183  break;
184  }
185  }
186 
187  if (channelSwitch == false) {
188  // ###############################################
189  // ### Induction Planes ###
190  // ###############################################
191 
192  //THE INDUCTION PLANE METHOD HAS NOT YET BEEN MODIFIED AND TESTED FOR REAL DATA.
193  // Or for detectors without a grid plane
194  //
195  if (sigType == geo::kInduction && !fSkipInd) {
196  threshold = fMinSigInd;
197  float negthr = -1.0 * threshold;
198  unsigned int bin = 1;
199  float minadc = 0;
200 
201  // find the dips
202  while (bin < (fDataSize - 1)) { // loop over ticks
203  float thisadc = holder[bin];
204  float nextadc = holder[bin + 1];
205  if (thisadc < negthr &&
206  nextadc < negthr) { // new region, require two ticks above threshold
207  // step back to find zero crossing
208  unsigned int place = bin;
209  while (thisadc <= 0 && bin > 0) {
210  bin--;
211  thisadc = holder[bin];
212  }
213  float hittime = bin + thisadc / (thisadc - holder[bin + 1]);
214  maxTimes.push_back(hittime);
215 
216  // step back more to find the hit start time
217  uint32_t stop;
218  if (fIndCutoff < (int)bin) { stop = bin - fIndCutoff; }
219  else {
220  stop = 0;
221  }
222  while (thisadc < threshold && bin > stop) {
223  bin--;
224  thisadc = holder[bin];
225  }
226  if (bin >= 2) bin -= 2;
227  while (thisadc > threshold && bin > stop) {
228  bin--;
229  thisadc = holder[bin];
230  }
231  startTimes.push_back(bin + 1);
232  // now step forward from hit time to find end time, area of dip
233  bin = place;
234  thisadc = holder[bin];
235  minadc = thisadc;
236  bin++;
237  totSig = fabs(thisadc);
238  while (thisadc < negthr && bin < fDataSize) {
239  totSig += fabs(thisadc);
240  thisadc = holder[bin];
241  if (thisadc < minadc) minadc = thisadc;
242  bin++;
243  }
244  endTimes.push_back(bin - 1);
245  peakHeight.push_back(-1.0 * minadc);
246  charge.push_back(totSig);
247  hitrms.push_back(5.0);
248 
249  // don't look for a new hit until it returns to baseline
250  while (thisadc < 0 && bin < fDataSize) {
251  bin++;
252  if (bin == fDataSize) break;
253  thisadc = holder[bin];
254  }
255  } // end region
256  bin++;
257  } // loop over ticks
258  }
259 
260  // ###############################################
261  // ### Collection Plane ###
262  // ###############################################
263 
264  else if (sigType == geo::kCollection) {
265  threshold = fMinSigCol;
266 
267  float madc = threshold;
268  int ibin = 0;
269  int start = 0;
270  int end = 0;
271  unsigned int bin = 0;
272 
273  while (bin < fDataSize) {
274  float thisadc = holder[bin];
275  madc = threshold;
276  ibin = 0;
277 
278  if (thisadc > madc) {
279  start = bin;
280 
281  if (thisadc > threshold && bin < fDataSize) {
282  while (thisadc > threshold && bin < fDataSize) {
283  if (thisadc > madc) {
284  ibin = bin;
285  madc = thisadc;
286  }
287  bin++;
288  if (bin == fDataSize) break;
289  thisadc = holder[bin];
290  }
291  }
292  else {
293  bin++;
294  }
295 
296  end = bin - 1;
297 
298  if (start != end) {
299  maxTimes.push_back(ibin);
300  peakHeight.push_back(madc);
301  startTimes.push_back(start);
302  endTimes.push_back(end);
303 
304  totSig = 0;
305  myrms = 0;
306  mynorm = 0;
307 
308  int moreTail = std::ceil(fIncludeMoreTail * (end - start));
309  if (moreTail < fColMinWindow) moreTail = fColMinWindow;
310 
311  for (int i = start - moreTail; i <= end + moreTail; i++) {
312  if (i < (int)(holder.size()) && i >= 0) {
313  float temp = ibin - i;
314  myrms += temp * temp * holder[i];
315 
316  totSig += holder[i];
317  }
318  }
319 
320  charge.push_back(totSig);
321  mynorm = totSig;
322  myrms /= mynorm;
323  hitrms.push_back(sqrt(myrms));
324  }
325  }
326  start = 0;
327  end = 0;
328  bin++;
329  }
330  }
331  }
332 
333  int numHits(0); //NUMBER OF CONSECUTIVE HITS BEING FITTED.
334  int hitIndex(0); //INDEX OF CURRENT HIT IN SEQUENCE.
335  double amplitude(0), position(0); //FIT PARAMETERS.
336  double start(0), end(0);
337  double amplitudeErr(0), positionErr(0); //FIT ERRORS.
338  double goodnessOfFit(0), chargeErr(0); //CHI2/NDF and error on charge.
339  double hrms(0);
340 
341  numHits = maxTimes.size();
342  for (int i = 0; i < numHits; ++i) {
343  amplitude = peakHeight[i];
344  position = maxTimes[i];
345  start = startTimes[i];
346  end = endTimes[i];
347  hrms = hitrms[i];
348  amplitudeErr = -1;
349  positionErr = 1.0;
350  goodnessOfFit = -1;
351  chargeErr = -1;
352  totSig = charge[i];
353 
354  std::vector<geo::WireID> wids = wireReadoutGeom.ChannelToWire(channel);
355  geo::WireID wid = wids[0];
356 
357  if (start >= end) {
358  mf::LogWarning("RawHitFinder_module")
359  << "Hit start " << start << " is >= hit end " << end;
360  continue;
361  }
362 
363  recob::HitCreator hit(*digitVec, //RAW DIGIT REFERENCE.
364  wid, //WIRE ID.
365  start, //START TICK.
366  end, //END TICK.
367  hrms, //RMS.
368  position, //PEAK_TIME.
369  positionErr, //SIGMA_PEAK_TIME.
370  amplitude, //PEAK_AMPLITUDE.
371  amplitudeErr, //SIGMA_PEAK_AMPLITUDE.
372  totSig, //HIT_INTEGRAL.
373  chargeErr, //HIT_SIGMA_INTEGRAL.
374  std::accumulate(holder.begin() + (int)start,
375  holder.begin() + (int)end,
376  0.), //SUMMED CHARGE.
377  std::accumulate(holder.begin() + (int)start,
378  holder.begin() + (int)end,
379  0.), //SUMMED CHARGE. TO BE MODIFIED
380  1, //MULTIPLICITY.
381  -1, //LOCAL_INDEX.
382  goodnessOfFit, //WIRE ID.
383  int(end - start + 1) //DEGREES OF FREEDOM.
384  );
385  hcol.emplace_back(hit.move(), digitVec);
386 
387  ++hitIndex;
388  }
389  }
390 
391  hcol.put_into(evt);
392  }
393 
395 
396 } // end namespace hit
float GetPedestal() const
Definition: RawDigit.h:221
const ADCvector_t & ADCs() const
Reference to the compressed ADC count vector.
Definition: RawDigit.h:209
void produce(art::Event &evt) override
ULong64_t Samples() const
Number of samples in the uncompressed ADC data.
Definition: RawDigit.h:217
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
ChannelID_t Channel() const
DAQ channel this raw data was read from.
Definition: RawDigit.h:213
std::vector< double > fAreaNorms
Definition of basic raw digits.
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
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:261
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
Class managing the creation of a new recob::Hit object.
Definition: HitCreator.h:87
Helper functions to create a hit.
A class handling a collection of hits and its associations.
Definition: HitCreator.h:489
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
Signal from induction planes.
Definition: geo_types.h:147
void reserve(size_t new_size)
Prepares the collection to host at least new_size hits.
Definition: HitCreator.h:608
enum geo::_plane_sigtype SigType_t
Enumerate the possible plane projections.
Collect all the RawData header files together.
RawHitFinder(fhicl::ParameterSet const &pset)
void emplace_back(recob::Hit &&hit, art::Ptr< recob::Wire > const &wire=art::Ptr< recob::Wire >(), art::Ptr< raw::RawDigit > const &digits=art::Ptr< raw::RawDigit >())
Adds the specified hit to the data collection.
Definition: HitCreator.cxx:299
float bin[41]
Definition: plottest35.C:14
Definition of data types for geometry description.
void put_into(art::Event &)
Moves the data into an event.
Definition: HitCreator.h:622
Detector simulation of raw signals on wires.
ProducesCollector & producesCollector() noexcept
raw::Compress_t Compression() const
Compression algorithm used to store the ADC counts.
Definition: RawDigit.h:229
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::string fCalDataModuleLabel
art::InputTag fDigitModuleLabel
TCEvent evt
Definition: DataStructs.cxx:8
void Uncompress(const std::vector< short > &adc, std::vector< short > &uncompressed, raw::Compress_t compress)
Uncompresses a raw data buffer.
Definition: raw.cxx:744
Signal from collection planes.
Definition: geo_types.h:148