AdcNoiseSignalFinder_tool.cc
Go to the documentation of this file.
1 // AdcNoiseSignalFinder_tool.cc
2 
3 #include "AdcNoiseSignalFinder.h"
4 #include <iostream>
5 
6 using std::string;
7 using std::cout;
8 using std::endl;
9 
10 //**********************************************************************
11 // Local defintitions.
12 //**********************************************************************
13 
14 namespace {
15 
16 string boolToString(bool val) {
17  return val ? "true" : "false";
18 }
19 
20 } // end unnamed namespace
21 
22 //**********************************************************************
23 // Class methods.
24 //**********************************************************************
25 
27 : m_LogLevel(ps.get<int>("LogLevel")),
28  m_SigFracMax(ps.get<float>("SigFracMax")),
29  m_ThresholdMin(ps.get<float>("ThresholdMin")),
30  m_ThresholdRatio(ps.get<float>("ThresholdRatio")),
31  m_ThresholdRatioTol(ps.get<float>("ThresholdRatioTol")),
32  m_MaxLoop(ps.get<unsigned int>("MaxLoop")),
33  m_BinsBefore(ps.get<unsigned int>("BinsBefore")),
34  m_BinsAfter(ps.get<unsigned int>("BinsAfter")),
35  m_FlagPositive(ps.get<bool>("FlagPositive")),
36  m_FlagNegative(ps.get<bool>("FlagNegative")) {
37  const string myname = "AdcNoiseSignalFinder::ctor: ";
38  if ( m_LogLevel >= 1 ) {
39  cout << myname << "Configuration parameters:" << endl;
40  cout << myname << " LogLevel: " << m_LogLevel << endl;
41  cout << myname << " SigFracMax: " << m_SigFracMax << endl;
42  cout << myname << " ThresholdMin: " << m_ThresholdMin << endl;
43  cout << myname << " ThresholdRatio: " << m_ThresholdRatio << endl;
44  cout << myname << " ThresholdRatioTol: " << m_ThresholdRatioTol << endl;
45  cout << myname << " BinsBefore: " << m_BinsBefore << endl;
46  cout << myname << " MaxLoop: " << m_MaxLoop << endl;
47  cout << myname << " BinsAfter: " << m_BinsAfter << endl;
48  cout << myname << " FlagPositive: " << boolToString(m_FlagPositive) << endl;
49  cout << myname << " FlagNegative: " << boolToString(m_FlagNegative) << endl;
50  }
51 }
52 
53 //**********************************************************************
54 
56  const string myname = "AdcNoiseSignalFinder::update: ";
57  DataMap ret;
58  AdcIndex nsam = acd.samples.size();
59  if ( m_ThresholdRatio <= 0 ) {
60  cout << myname << "Invalid noise ratio: " << m_ThresholdRatio << endl;
61  return ret.setStatus(2);
62  }
63  if ( m_ThresholdRatioTol <= 0 ) {
64  cout << myname << "Invalid noise ratio tolerance: " << m_ThresholdRatioTol << endl;
65  return ret.setStatus(2);
66  }
67  if ( nsam == 0 ) {
68  cout << myname << "ERROR: No samples found in channel " << acd.channel() << endl;
69  acd.signal.clear();
70  acd.rois.clear();
71  return ret.setStatus(1);
72  }
73  if ( m_LogLevel >= 3 ) cout << myname << "Finding ROIs for channel " << acd.channel() << endl;
74  float thr = m_ThresholdMin;
75  float noise = 0.0;
76  float sigfrac = 0.0;
77  Index nloop = 0;
78  float trtgt = m_ThresholdRatio;
79  float trmin = trtgt - m_ThresholdRatioTol;
80  float trmax = trtgt + m_ThresholdRatioTol;
81  float thrTooLow = m_ThresholdMin; // We know threshold is at or above this value
82  float thrTooHigh = 1.e20; // We know threshold is below this value
83  while ( true ) {
84  // Find the and flag the signal.
85  AdcIndex nsamlo = m_BinsBefore;
86  AdcIndex nsamhi = m_BinsAfter;
87  AdcIndex nsamAbove = 0;
88  AdcIndex isamUnknown = 0; // First sample not known to be in or outside a ROI.
89  acd.signal.clear();
90  acd.signal.resize(nsam, false);
91  for ( AdcIndex isam=0; isam<nsam; ++isam ) {
92  bool keep = ( m_FlagPositive && acd.samples[isam] > thr ) ||
93  ( m_FlagNegative && acd.samples[isam] < -thr );
94  if ( keep ) {
95  ++nsamAbove;
96  AdcIndex jsam1 = isam > nsamlo ? isam - nsamlo : 0;
97  if ( jsam1 < isamUnknown ) jsam1 = isamUnknown;
98  AdcIndex jsam2 = isam + nsamhi + 1;
99  if ( jsam2 > nsam ) jsam2 = nsam;
100  if ( m_LogLevel >= 5 ) cout << myname << "Trigger: " << isam << ", range: ["
101  << jsam1 << ", " << jsam2 << ")" << endl;
102  for ( AdcIndex jsam=jsam1; jsam<jsam2; ++jsam ) acd.signal[jsam] = true;
103  isamUnknown = jsam2;
104  }
105  }
106  // Evaluate the noise and signal fraction.
107  Index nsig = 0;
108  Index nnsg = 0;
109  float ssqsum = 0.0;
110  for ( AdcIndex isam=0; isam<nsam; ++isam ) {
111  if ( acd.signal[isam] ) {
112  ++nsig;
113  } else {
114  ++nnsg;
115  float val = acd.samples[isam];
116  ssqsum += val*val;
117  }
118  }
119  noise = nnsg ? sqrt(ssqsum/nnsg) : 0.0;
120  sigfrac = float(nsig)/nsam;
121  // Check is we need andother loop and rest threshold accordingly.
122  ++nloop;
123  if ( m_LogLevel >= 4 ) {
124  cout << myname << " Loop " << nloop << endl;
125  cout << myname << " Threshold: " << thr << endl;
126  cout << myname << " Noise: " << noise << endl;
127  cout << myname << " Thr/noise: " << (noise > 0.0 ? thr/noise : -999) << endl;
128  cout << myname << " # ticks above threshold: " << nsamAbove << endl;
129  cout << myname << " SigFrac: " << sigfrac << endl;
130  acd.roisFromSignal();
131  cout << myname << " # ROI: " << acd.rois.size() << endl;
132  }
133  if ( nloop >= m_MaxLoop ) {
134  cout << myname << "WARNING: Channel " << acd.channel() << " exiting after "
135  << nloop << " loops." << endl;
136  break;
137  }
138  // Use current noise estimate to evaluate the target threshold and range.
139  float thrtgt = trtgt*noise;
140  float thrmin = trmin*noise;
141  float thrmax = trmax*noise;
142  float thrOld = thr;
143  // Noise level too high ==> raise threshold.
144  if ( sigfrac > m_SigFracMax || thrmin > thr ) {
145  thrTooLow = thr;
146  // Double the threshold unless that goes past the known upper limit.
147  // In that case, go halfway.
148  float thrEst = thrtgt > thr ? 2.0*thrtgt : 2.0*thr;
149  if ( thrEst < thrTooHigh ) thr = thrEst;
150  else thr = 0.5*(thr + thrTooHigh);
151  } else if ( thr <= m_ThresholdMin ) {
152  break;
153  } else if ( thrmax >= thr ) {
154  break;
155  } else {
156  // Noise level too low ==> lower threshold.
157  thrTooHigh = thr;
158  // Use the target as the next estimate unless it below the thw known lower
159  // limit. In the latter case, go halfway to that limit.
160  if ( thrtgt > thrTooLow ) thr = thrtgt;
161  else thr = 0.5*(thr + thrTooLow);
162  }
163  float dthrmax = 0.001*thrOld;
164  if ( m_LogLevel >= 1 ) {
165  if ( fabs(thr - thrOld) < dthrmax ) {
166  Index chstat = acd.channelStatus();
167  if ( m_LogLevel >= 2 || (m_LogLevel == 1 && chstat == 0) ) {
168  string schstat = chstat==0 ? "Good" :
169  chstat==1 ? "Bad" :
170  chstat==2 ? "Noisy" : "UnkownStatus";
171  cout << myname << "WARNING: " << schstat << " channel " << acd.channel()
172  << " exiting prematurely due to threshold convergence." << endl;
173  }
174  break;
175  }
176  }
177  }
178  acd.roisFromSignal();
179  ret.setFloat( "nsfSigFrac", sigfrac);
180  ret.setFloat( "nsfNoise", noise);
181  ret.setFloat("nsfThreshold", thr);
182  ret.setInt( "nsfLoopCount", nloop);
183  ret.setInt( "nsfRoiCount", acd.rois.size());
184  if ( m_LogLevel >= 4 ) ret.print(myname);
185  acd.metadata[ "nsfSigFrac"] = sigfrac;
186  acd.metadata[ "nsfNoise"] = noise;
187  acd.metadata["nsfThreshold"] = thr;
188  acd.metadata["nsfLoopCount"] = nloop;
189  acd.metadata[ "nsfRoiCount"] = acd.rois.size();
190 
191  return ret;
192 }
193 //**********************************************************************
194 
196  AdcChannelData acdtmp;
197  acdtmp.samples = acd.samples;
198  return update(acdtmp);
199 }
200 
201 //**********************************************************************
202 
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
void setFloat(Name name, float val)
Definition: DataMap.h:133
G4double thr[100]
Definition: G4S2Light.cc:59
unsigned int Index
DataMap & setStatus(int stat)
Definition: DataMap.h:130
std::string string
Definition: nybbler.cc:12
void print(std::ostream *pout) const
Definition: DataMap.h:245
AdcRoiVector rois
AdcNoiseSignalFinder(fhicl::ParameterSet const &ps)
DataMap view(const AdcChannelData &acd) const override
void setInt(Name name, int val)
Definition: DataMap.h:131
static constexpr double ps
Definition: Units.h:99
unsigned int AdcIndex
Definition: AdcTypes.h:15
Index channelStatus() const
Channel channel() const
AdcFilterVector signal
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
int bool
Definition: qglobal.h:345
DataMap update(AdcChannelData &acd) const override
AdcSignalVector samples
QTextStream & endl(QTextStream &s)