Merge pull request #61 from sjgardiner/phaseI-beam-cuts

Phase I beam cuts

Author: brichards64

DataModel/ANNIEalgorithms.h

@@ -25,19 +25,19 @@ template<typename ElementType> void ComputeMeanAndVariance(
   }
 
   size_t num_samples = 0;
-  double mean_x2 = 0.;
+  double m2 = 0.;
   mean = 0.;
 
   for (const ElementType& x : data) {
     ++num_samples;
     double delta = x - mean;
     mean += delta / num_samples;
     double delta2 = x - mean;
-    mean_x2 = delta * delta2;
+    m2 += delta * delta2;
     if (num_samples == sample_cutoff) break;
   }
 
-  var = mean_x2 / (num_samples - 1);
+  var = m2 / (num_samples - 1);
   return;
 }
 

DataModel/BeamStatus.cpp

@@ -17,4 +17,18 @@ void BeamStatus::clear()
   time_ = TimeClass(0ull);
   pot_ = 0.;
   condition_ = BeamCondition::Missing;
+  data_.clear();
+  cuts_.clear();
+}
+
+void BeamStatus::add_measurement(const std::string& device_name,
+  uint64_t ms_since_epoch, const BeamDataPoint& bdp)
+{
+  data_[device_name] = std::pair<uint64_t, BeamDataPoint>(ms_since_epoch, bdp);
+}
+
+void BeamStatus::add_measurement(const std::string& device_name,
+  uint64_t ms_since_epoch, double value, const std::string& unit)
+{
+  add_measurement(device_name, ms_since_epoch, BeamDataPoint(value, unit));
 }

DataModel/BeamStatus.h

@@ -8,9 +8,11 @@
 // standard library includes
 #include <iostream>
 #include <fstream>
+#include <map>
 #include <sstream>
 
 // ToolAnalysis includes
+#include "BeamDataPoint.h"
 #include "SerialisableObject.h"
 #include "TimeClass.h"
 
@@ -59,6 +61,9 @@ class BeamStatus : public SerialisableObject {
     inline TimeClass time() const { return time_; }
     inline double pot() const { return pot_; }
     inline BeamCondition condition() const { return condition_; }
+    inline const std::map<std::string, std::pair<uint64_t, BeamDataPoint> >&
+      data() const { return data_; }
+    inline const std::map<std::string, bool>& cuts() const { return cuts_; }
 
     // Good beam minibuffers for the analysis should use is_beam() == true
     // and ok() == true. Good non-beam minibuffers for the analysis
@@ -73,14 +78,52 @@ class BeamStatus : public SerialisableObject {
       { return condition_ != BeamCondition::Bad
           && condition_ != BeamCondition::Missing; }
 
+    // Beam quality cut checks
+    inline bool passed_cut(const std::string& cut_name) const {
+      return cuts_.at(cut_name);
+    }
+    inline bool passed_all_cuts() const {
+      // Returns true if all of the entries in the cuts_ map have true
+      // boolean values, or false otherwise.
+      return std::all_of(cuts_.cbegin(), cuts_.cend(),
+        [](const std::pair<std::string, bool>& p) -> bool { return p.second; });
+    }
+
     inline void set_time(TimeClass time) { time_ = time; }
     inline void set_pot(double POT) { pot_ = POT; }
     inline void set_condition(BeamCondition bc) { condition_ = bc; }
 
+    // Adds a new monitoring device measurement to the data map
+    void add_measurement(const std::string& device_name,
+      uint64_t ms_since_epoch, const BeamDataPoint& bdp);
+
+    void add_measurement(const std::string& device_name,
+      uint64_t ms_since_epoch, double value, const std::string& unit);
+
+    // Adds a new cut result to the cuts map (passed == true means that the
+    // beam spill survived the quality cut)
+    // TODO: Note that this will overwrite an existing cut with the same
+    // name. Think more carefully about whether this is the optimal behavior.
+    inline void add_cut(const std::string& cut_name, bool passed) {
+      cuts_[cut_name] = passed;
+    }
+
     inline bool Print() {
       std::cout << "Timestamp : " << time_ << '\n';
       std::cout << "Beam Intensity [ETOR875] : " << pot_ << " POT\n";
       std::cout << "Condition : " << condition_ << '\n';
+      for (auto& outer_pair : data_) {
+        std::cout << "Monitoring Data for device : " << outer_pair.first
+          << '\n';
+        std::cout << "Measurement timestamp : " << outer_pair.second.first
+          << " ms since epoch\n";
+        outer_pair.second.second.Print();
+      }
+      for (auto& pair : cuts_) {
+        if (pair.second) std::cout << "PASSED ";
+        else std::cout << "FAILED ";
+        std::cout << " the " << pair.first << " beam quality cut\n";
+      }
       return true;
     }
 
@@ -106,18 +149,54 @@ class BeamStatus : public SerialisableObject {
     /// successfully, but the current beam spill should be ignored in the
     /// analysis). Reasons for marking a beam minibuffer as "Bad" include
     /// a very low POT value (suggesting that the beam monitor E:TOR875
-    /// was off at the time, Hefty mode information suggesting that one or
-    /// more self-trigger minibuffers might be missing after the current beam
-    /// spill, etc.
+    /// was off at the time), a low peak horn current, etc.
     BeamCondition condition_;
 
+    // Map storing data from the beam monitoring instruments for the spill
+    // of interest. Keys are device names, values are pairs where the
+    // first element is a timestamp for the measurement (in ms since the Unix
+    // epoch) and the second element is a BeamDataPoint object.
+    //
+    // ** Known device names (from the Intensity Frontier beam database) **
+    // E:1DCNT = Accumulating count of the number of "$1D" triggers (a "$1D"
+    //           signal indicates that the Booster is preparing to send a beam
+    //           pulse to ANNIE).
+    //
+    // E:HI860, E:VI860, E:HP860, E:VP860 = beam intensity (I) or position (P)
+    //                                      monitor reading at position 860.
+    //                                      Horizontal (H) or vertical (V)
+    //                                      position monitors are read
+    //                                      separately.
+    //
+    // E:HP875, E:VP875 = similar monitors at position 875
+    //
+    // E:HITG1, E:VITG1, E:VPTG1, E:HPTG1 = similar monitors at position TG1
+    //
+    // E:HPTG2, E:HITG2, E:VITG2, E:VPTG2 = similar monitors at position TG2
+    //
+    // E:THCURR = Horn current
+    // E:TOR860 = POT as measured by toroid at position 860
+    //            (further away from target than 875)
+    //
+    // E:TOR875 = POT as measured by toroid at position 875
+    std::map<std::string, std::pair<uint64_t, BeamDataPoint> > data_;
+
+    // Map storing the results of each of the beam quality cut checks.
+    // Keys are cut names, values are boolean variables telling whether
+    // the spill of interest passed (true) or failed (false) the beam
+    // quality cut. For good beam spills to be used for analysis, all of
+    // the entries in this map should have a true boolean value.
+    std::map<std::string, bool> cuts_;
+
     template<class Archive> void serialize(Archive & ar,
       const unsigned int version)
     {
       if (!serialise) return;
       ar & time_;
       ar & pot_;
       ar & condition_;
+      ar & data_;
+      ar & cuts_;
     }
 };
 #endif

UserTools/ADCCalibrator/ADCCalibrator.cpp

@@ -1,3 +1,5 @@
+#include <string>
+
 // ToolAnalysis includes
 #include "ADCCalibrator.h"
 #include "ANNIEalgorithms.h"
@@ -56,6 +58,9 @@ bool ADCCalibrator::Execute() {
     const auto& channel_key = temp_pair.first;
     const auto& raw_waveforms = temp_pair.second;
 
+    Log("Making calibrated waveforms for ADC channel " +
+      std::to_string(channel_key.GetDetectorElementIndex()), 3, verbosity);
+
     calibrated_waveform_map[channel_key] = make_calibrated_waveforms(
       raw_waveforms);
   }
@@ -74,17 +79,22 @@ void ADCCalibrator::ze3ra_baseline(
   const std::vector< Waveform<unsigned short> >& raw_data,
   double& baseline, double& sigma_baseline, size_t num_baseline_samples)
 {
-  // All F-distribution probabilities below this value will pass the
-  // variance consistency test in ze3ra_baseline()
-  double q_critical;
-  m_variables.Get("QCritical", q_critical);
+  int verbosity;
+  m_variables.Get("verbose", verbosity);
+
+  // All F-distribution probabilities above this value will pass the
+  // variance consistency test in ze3ra_baseline(). That is, p_critical
+  // is the maximum p-value for which we will reject the null hypothesis
+  // of equal variances.
+  double p_critical;
+  m_variables.Get("PCritical", p_critical);
 
   // Signal ADC means, variances, and F-distribution probability values
-  // ("Q") for the first num_baseline_samples from each minibuffer
+  // ("P") for the first num_baseline_samples from each minibuffer
   // (in Hefty mode) or from each sub-minibuffer (in non-Hefty mode)
   std::vector<double> means;
   std::vector<double> variances;
-  std::vector<double> Qs;
+  std::vector<double> Ps;
 
   // Hefty mode uses multiple minibuffers, non-Hefty mode uses a single
   // minibuffer
@@ -136,11 +146,20 @@ void ADCCalibrator::ze3ra_baseline(
     else F = sigma2_jp1 / sigma2_j;
 
     double nu = (num_baseline_samples - 1) / 2.;
-    double Q = std::tgamma(2*nu)
-      * annie_math::Incomplete_Beta_Function(1. / (1. + F), nu, nu)
-      / (2. * std::tgamma(nu));
+    double P = annie_math::Regularized_Beta_Function(1. / (1. + F), nu, nu);
+
+    // Two-tailed hypothesis test (we need to exclude unusually small values
+    // as well as unusually large ones). The tails have equal sizes, so we
+    // may use symmetry and simply multiply our earlier result by 2.
+    P *= 2.;
+
+    // I've never seen this problem (the numerical values for the regularized
+    // beta function that I've checked all fall within [0,1]), but the book
+    // Numerical Recipes includes this check in a similar block of code,
+    // so I'll add it just in case.
+    if (P > 1.) P = 2. - P;
 
-    Qs.push_back(Q);
+    Ps.push_back(P);
   }
 
   // Compute the mean and standard deviation of the baseline signal
@@ -149,30 +168,61 @@ void ADCCalibrator::ze3ra_baseline(
   // the critical value.
   baseline = 0.;
   sigma_baseline = 0.;
+  double variance_baseline = 0.;
   size_t num_passing = 0;
-  for (size_t k = 0; k < Qs.size(); ++k) {
-    if (Qs.at(k) < q_critical) {
+  for (size_t k = 0; k < Ps.size(); ++k) {
+    if (Ps.at(k) > p_critical) {
       ++num_passing;
       baseline += means.at(k);
-      sigma_baseline += std::sqrt( variances.at(k) );
+      variance_baseline += variances.at(k);
     }
   }
 
-  if (num_passing > 0) {
+  if (num_passing > 1) {
     baseline /= num_passing;
-    sigma_baseline /= num_passing;
+
+    variance_baseline *= static_cast<double>(num_baseline_samples - 1)
+      / (num_passing*num_baseline_samples - 1);
+    // Now that we've combined the sample variances correctly, take the
+    // square root to get the standard deviation
+    sigma_baseline = std::sqrt( variance_baseline );
+  }
+  else if (num_passing == 1) {
+    // We only have one set of sample statistics, so all we need to
+    // do is take the square root of the variance to get the standard
+    // deviation.
+    sigma_baseline = std::sqrt( variance_baseline );
   }
   else {
     // If none of the minibuffers passed the F-distribution test,
-    // choose the one closest to passing and adopt its baseline statistics
+    // choose the one closest to passing (i.e., the one with the largest
+    // P-value) and adopt its baseline statistics. For a sufficiently large
+    // number of minibuffers (e.g., 40), such a situation should be very rare.
     // TODO: consider changing this approach
-    auto min_iter = std::min_element(Qs.cbegin(), Qs.cend());
-    int min_index = std::distance(Qs.cbegin(), min_iter);
+    auto max_iter = std::max_element(Ps.cbegin(), Ps.cend());
+    int max_index = std::distance(Ps.cbegin(), max_iter);
 
-    baseline = means.at(min_index);
-    sigma_baseline = std::sqrt( variances.at(min_index) );
+    baseline = means.at(max_index);
+    sigma_baseline = std::sqrt( variances.at(max_index) );
   }
 
+  std::string mb_temp_string = "minibuffer";
+  if ( !hefty_mode ) mb_temp_string = "sub-minibuffer";
+
+  if (verbosity >= 4) {
+    for ( size_t x = 0; x < Ps.size(); ++x ) {
+      Log("  " + mb_temp_string + " " + std::to_string(x) + ", mean = "
+        + std::to_string(means.at(x)) + ", var = "
+        + std::to_string(variances.at(x)) + ", p-value = "
+        + std::to_string(Ps.at(x)), 4, verbosity);
+    }
+  }
+
+  Log(std::to_string(num_passing) + " " + mb_temp_string + " pairs passed the"
+    " F-test", 3, verbosity);
+  Log("Baseline estimate: " + std::to_string(baseline) + " ± "
+    + std::to_string(sigma_baseline) + " ADC counts", 3, verbosity);
+
 }
 
 std::vector< CalibratedADCWaveform<double> >

UserTools/ADCCalibrator/ADCCalibrator.h

@@ -25,7 +25,7 @@ class ADCCalibrator : public Tool {
     /// @brief Compute the baseline for a particular RawChannel
     /// object using a technique taken from the ZE3RA code.
     /// @details See section 2.2 of https://arxiv.org/pdf/1106.0808.pdf for a
-    /// full description of the algorithm.
+    /// description of the algorithm.
     void ze3ra_baseline(const std::vector< Waveform<unsigned short> >& raw_data,
       double& baseline, double& sigma_baseline, size_t num_baseline_samples);