HighEnergyRemover.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file: HighEnergyRemover.cxx
3 // \brief Remove hits associated with high-energy activity.
4 // \author Justin Vasel <justin.vasel@gmail.com>
5 // \date 2020-08-09
6 ////////////////////////////////////////////////////////////////////////
7 
8 // STL includes
9 #include <math.h>
10 
11 // ROOT includes
12 #include "TF1.h"
13 #include "TH1.h"
14 #include "TMath.h"
15 #include "TSpectrum.h"
16 
17 // ART includes
19 
20 // NOvASoft includes
21 #include "RecoBase/CellHit.h"
24 
25 // ............................................................................
27 fHitVetoMap(&vetoMapHit),
28 fNumRemoved(0),
29 fPeakThreshold(5),
30 fVetoWindowPadding(0.25)
31 {
32 }
33 
34 
35 // ............................................................................
37 
38 
39 // ............................................................................
41 {
42  int32_t sum = 0;
43  for (std::pair<int32_t, int32_t> w : this->fTDCVetoWindows) {
44  sum += (w.second - w.first);
45  }
46  return sum;
47 }
48 
49 
50 // ............................................................................
51 std::pair<int, int> sn::HighEnergyRemover::RegionAroundPeak(TH1F* h, int peakBin, double baseline, float padding)
52 {
53  int binLow = peakBin;
54  int binHigh = peakBin;
55 
56  // Walk lower until we drop below baseline
57  while (true) {
58  if (binLow==1) break;
59 
60  binLow -= 1;
61  double y = h->GetBinContent(binLow);
62 
63  if (y<baseline) break;
64  }
65 
66  // Walk higher until we drop below baseline
67  while (true) {
68  if (binHigh==h->GetNbinsX()) break;
69 
70  binHigh += 1;
71  double y = h->GetBinContent(binHigh);
72 
73  if (y<baseline) break;
74  }
75 
76  // Get bin-span between low and high bins and add padding
77  int binSpan = binHigh - binLow;
78  int binPadding = std::ceil(padding * binSpan);
79  int32_t xLo = h->GetBinCenter(binLow-binPadding);
80  int32_t xHi = h->GetBinCenter(binHigh+binPadding);
81 
82  return std::make_pair<int32_t, int32_t>((int32_t)xLo, (int32_t)xHi);
83 }
84 
85 
86 // ............................................................................
88 {
89  for (std::pair<int32_t, int32_t> w : this->fTDCVetoWindows) {
90  if (h.TDC() >= w.first && h.TDC() <= w.second) return true;
91  }
92  return false;
93 }
94 
95 
96 // ............................................................................
98 {
99  return this->fTDCVetoWindows.size();
100 }
101 
102 
103 // ............................................................................
105 {
106  int sum = 0;
107  for (std::pair<int32_t, int32_t> w : this->fTDCVetoWindows) {
108  sum += w.second - w.first;
109  }
110  return sum;
111 }
112 
113 
114 // ............................................................................
116 {
117  return this->fTotalTDC;
118 }
119 
120 
121 // ............................................................................
123 {
124  /* Compute binning */
125  // NOTE: 64 TDC = 1 us
126  int32_t binWidth = 64;
127 
128  int32_t tdc0 = hits.front()->TDC();
129  int32_t tdc1 = hits.back()->TDC();
130  this->fTotalTDC = tdc1 - tdc0;
131 
132  int32_t nBins = std::ceil((float)this->fTotalTDC / (float)binWidth);
133  int32_t bin0 = tdc0;
134  int32_t bin1 = tdc0 + nBins * binWidth;
135 
136  /* Fill histogram */
137  TH1F *h = new TH1F("h", "", nBins, bin0, bin1);
138  for (art::Ptr<rb::CellHit> aHit : hits) {
139  if (this->fHitVetoMap->ContainsHit(*aHit)) continue;
140  h->Fill(aHit->TDC());
141  }
142 
143  /* Search for peaks with sigma=2 and threshold=0.5 */
144  int maxPeaks = 100;
145  int nPeaksFound = 0;
146  TSpectrum *s = new TSpectrum(maxPeaks);
147  try {
148  nPeaksFound = s->Search(h, 2, "goff", 0.5);
149  } catch (...) {
150  std::cout << "• Warning: Max number of peaks (" << maxPeaks << ") exceeded. Will not look for any more." << std::endl;
151  }
152  double *xpeaks = s->GetPositionX();
153 
154  /* Fit a constant (0-degree polynomial) to the histogram */
155  // NOTE: don't include bins at edges as they can make the fit worse
156  TF1 *fline = new TF1("fline", "pol0", bin0+2, bin1-2);
157  h->Fit("fline", "Lnq");
158 
159  /* Loop over peaks and determine which should be kept */
160  int nPeaksPassed = 0;
161  for (int pIdx=0; pIdx<nPeaksFound; pIdx++) {
162  double xp = xpeaks[pIdx];
163  int bin = h->GetXaxis()->FindBin(xp);
164  float baseline = fline->Eval(xp);
165  double yp = h->GetBinContent(bin);
166 
167  // If the height of this peak isn't large enough, we won't keep it
168  if (yp < baseline + this->fPeakThreshold * TMath::Sqrt(baseline)) continue;
169 
170  nPeaksPassed++;
171  std::cout << "• Peak " << nPeaksPassed << " at TDC=" << xp << ", y=" << yp << ", baseline=" << baseline << std::endl;
172 
173  // Find the low and high bins that encompass the region surrounding this peak
174  // with a user-defined padding on either side and add to TDC veto list
175  this->fTDCVetoWindows.push_back(this->RegionAroundPeak(h, bin, fline->Eval(xp), this->fVetoWindowPadding));
176  }
177 
178  /* Add cell hits from the veto regions to veto map */
179  for (art::Ptr<rb::CellHit> aHit : hits) {
180  if (this->fHitVetoMap->ContainsHit(*aHit)) continue;
181  if (!(this->HitInTDCVetoWindow(*aHit))) continue;
182 
183  this->fHitVetoMap->AddHit(*aHit);
184  }
185 
186  delete h;
187  delete s;
188  delete fline;
189 
190  return;
191 }
int nBins
Definition: plotROC.py:16
std::vector< std::pair< int32_t, int32_t > > fTDCVetoWindows
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
HighEnergyRemover(sn::HitVetoMap &vetoMapHit)
bool ContainsHit(rawdata::RawDigit h)
Definition: HitVetoMap.cxx:41
const XML_Char * s
Definition: expat.h:262
void hits()
Definition: readHits.C:15
void AddHit(rawdata::RawDigit h)
Definition: HitVetoMap.cxx:27
std::pair< int, int > RegionAroundPeak(TH1F *h, int peakBin, double baseline, float padding)
float bin[41]
Definition: plottest35.C:14
void remove(std::vector< art::Ptr< rb::CellHit >> &hits)
OStream cout
Definition: OStream.cxx:6
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
sn::HitVetoMap * fHitVetoMap
bool HitInTDCVetoWindow(rb::CellHit h)
Double_t sum
Definition: plot.C:31
Float_t w
Definition: plot.C:20
fvar< T > ceil(const fvar< T > &x)
Definition: ceil.hpp:11