LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
TruncMean.cxx
Go to the documentation of this file.
1 #ifndef TRUNCMEAN_CXX
2 #define TRUNCMEAN_CXX
3 
4 #include "TruncMean.h"
5 
6 
7 float TruncMean::CalcIterativeTruncMean(std::vector<float> v, const size_t& nmin,
8  const size_t& nmax, const size_t& currentiteration,
9  const size_t& lmin,
10  const float& convergencelimit,
11  const float& nsigma, const float& oldmed)
12 {
13 
14  auto const& mean = Mean(v);
15  auto const& med = Median(v);
16  auto const& rms = RMS(v);
17 
18  // if the vector length is below the lower limit -> return
19  if (v.size() < lmin)
20  return mean;
21 
22  // if we have passed the maximum number of iterations -> return
23  if (currentiteration >= nmax)
24  return mean;
25 
26  // if we passed the minimum number of iterations and the mean is close enough to the old value
27  float fracdiff = fabs(med-oldmed) / oldmed;
28  if ( (currentiteration >= nmin) && (fracdiff < convergencelimit) )
29  return mean;
30 
31  // if reached here it means we have to go on for another iteration
32 
33  // cutoff tails of distribution surrounding the mean
34  // use erase-remove : https://en.wikipedia.org/wiki/Erase%E2%80%93remove_idiom
35  // https://stackoverflow.com/questions/17270837/stdvector-removing-elements-which-fulfill-some-conditions
36  v.erase( std::remove_if( v.begin(), v.end(),
37  [med,nsigma,rms](const float& x) { return ( (x < (med-nsigma*rms)) || (x > (med+nsigma*rms)) ); }), // lamdda condition for events to be removed
38  v.end());
39 
40  return CalcIterativeTruncMean(v, nmin, nmax, lmin, currentiteration+1, convergencelimit, nsigma, med);
41 }
42 
43 void TruncMean::CalcTruncMeanProfile(const std::vector<float>& rr_v, const std::vector<float>& dq_v,
44  std::vector<float>& dq_trunc_v, const float& nsigma)
45 {
46 
47  // how many points to sample
48  int Nneighbor = (int)(_rad * 3 * 2);
49 
50  dq_trunc_v.clear();
51  dq_trunc_v.reserve( rr_v.size() );
52 
53  int Nmax = dq_v.size()-1;
54 
55  for (size_t n=0; n < dq_v.size(); n++) {
56 
57  // current residual range
58  float rr = rr_v.at(n);
59 
60  int nmin = n - Nneighbor;
61  int nmax = n + Nneighbor;
62 
63  if (nmin < 0) nmin = 0;
64  if (nmax > Nmax) nmax = Nmax;
65 
66  // vector for local dq values
67  std::vector<float> dq_local_v;
68 
69  for (int i=nmin; i < nmax; i++) {
70 
71  float dr = rr - rr_v[i];
72  if (dr < 0) dr *= -1;
73 
74  if (dr > _rad) continue;
75 
76  dq_local_v.push_back( dq_v[i] );
77 
78  }// for all ticks we want to scan
79 
80  if (dq_local_v.size() == 0) {
81  dq_trunc_v.push_back( dq_v.at(n) );
82  continue;
83  }
84 
85  // calculate median and rms
86  float median = Median(dq_local_v);
87  float rms = RMS(dq_local_v);
88 
89  float truncated_dq = 0.;
90  int npts = 0;
91  for (auto const& dq : dq_local_v) {
92  if ( ( dq < (median+rms * nsigma) ) && ( dq > (median-rms * nsigma) ) ){
93  truncated_dq += dq;
94  npts += 1;
95  }
96  }
97 
98  dq_trunc_v.push_back( truncated_dq / npts );
99  }// for all values
100 
101  return;
102 }
103 
104 float TruncMean::Mean(const std::vector<float>& v)
105 {
106 
107  float mean = 0.;
108  for (auto const& n : v) mean += n;
109  mean /= v.size();
110 
111  return mean;
112 }
113 
114 float TruncMean::Median(const std::vector<float>& v)
115 {
116 
117  if (v.size() == 1) return v[0];
118 
119  std::vector<float> vcpy = v;
120 
121  std::sort(vcpy.begin(), vcpy.end());
122 
123  float median = vcpy[ vcpy.size() / 2 ];
124 
125  return median;
126 }
127 
128 float TruncMean::RMS(const std::vector<float>& v)
129 {
130 
131  float avg = 0.;
132  for (auto const& val : v) avg += val;
133  avg /= v.size();
134  float rms = 0.;
135  for (auto const& val : v) rms += (val-avg)*(val-avg);
136  rms = sqrt( rms / ( v.size() - 1 ) );
137 
138  return rms;
139 }
140 
141 #endif
Float_t x
Definition: compare.C:6
Class def header for a class TruncMean.
float Median(const std::vector< float > &v)
Definition: TruncMean.cxx:114
double _rad
Definition: TruncMean.h:99
float RMS(const std::vector< float > &v)
Definition: TruncMean.cxx:128
double mean(const std::vector< short > &wf, size_t start, size_t nsample)
Definition: UtilFunc.cxx:15
float CalcIterativeTruncMean(std::vector< float > v, const size_t &nmin, const size_t &nmax, const size_t &currentiteration, const size_t &lmin, const float &convergencelimit, const float &nsigma, const float &oldmed=kINVALID_FLOAT)
Iteratively calculate the truncated mean of a distribution std::vector<float> v -> vector of values ...
Definition: TruncMean.cxx:7
float Mean(const std::vector< float > &v)
Definition: TruncMean.cxx:104
void CalcTruncMeanProfile(const std::vector< float > &rr_v, const std::vector< float > &dq_v, std::vector< float > &dq_trunc_v, const float &nsigma=1)
Given residual range and dq vectors return truncated local dq. Input vectors are assumed to be match ...
Definition: TruncMean.cxx:43
Char_t n[5]