Merge lp:~m-a-uchida/maus/GlobalDevelopment into lp:maus

[Melissa Uchida]

Just a very small change to track matching to improve the user experience, putting exceptions into debug mode.

=== modified file 'doc/doc_src/maus_user_guide.tex'
--- doc/doc_src/maus_user_guide.tex	2017-05-15 14:17:05 +0000
+++ doc/doc_src/maus_user_guide.tex	2017-06-08 17:10:31 +0000
@@ -60,7 +60,7 @@
 
 \input{reconstruction/globaltrackmatching/globaltrackmatching.tex}
 \input{reconstruction/globalpid/globalpid.tex}
-% \input{reconstruction/globalds/globalds.tex}
+\input{reconstruction/globalds/globalds.tex}
 
 \input{envelope_tool} % envelope utility - could be an appendix
 

=== added directory 'doc/doc_src/reconstruction/globalds'
=== added file 'doc/doc_src/reconstruction/globalds/globalds.tex'
--- doc/doc_src/reconstruction/globalds/globalds.tex	1970-01-01 00:00:00 +0000
+++ doc/doc_src/reconstruction/globalds/globalds.tex	2017-06-08 17:10:31 +0000
@@ -0,0 +1,91 @@
+\chapter{Accessing Global Tracks}
+\label{chapter:globalds}
+
+The global datastructure is set up so that the state of tracks at any point in the processing chain, import,
+track matching, pid, and track fitting, can easily be accessed in order to understand where a given track
+came from and what lead to its current state.
+
+To achieve this, tracks are packaged together in PrimaryChain objects which retain copies of the
+tracks at every step along the way. There will generally be one primary chain for every tracker track that
+has been matched to at least one additional detector, plus possible through primary chains for tracks
+that have been matched through the absorber.
+
+\section{The PrimaryChain Object}\label{sec:primarychain}
+
+A primary chain (\texttt{MAUS::DataStructure::Global::PrimaryChain}) has a number of member functions
+for accessing parameters describing the chain as well as the variety of tracks that are contained in it:
+
+\begin{description}
+  \item [get\_chain\_type()] A chain can hold either upstream tracks, downstream tracks, or through tracks. The chain type is an enumerator with the following values
+  \begin{description}
+    \item [kNoChainType] Used for initialization only
+    \item [kUS] A chain containing upstream tracks. This value furthermore implies that a track in the chain has been matched into a through chain.
+    \item [kDS] A chain containing downstream tracks. This value furthermore implies that a track in the chain has been matched into a through chain.
+    \item [kUSOrphan] A chain containing upstream tracks none of which has been matched into a through track
+    \item [kDSOrphan] A chain containing downstream tracks none of which has been matched into a through track
+    \item [kThrough] A chain containing through tracks. It will also reference daughter upstream and downstream chains
+  \end{description}
+  \item [get\_chain\_multiplicity] This is only a meaningful property for through chains. It describes whether this and another through chain are mutually exclusive
+        because they share daughter chains, for example if the same upstream track has been matched to multiple downstream tracks. The possible values of the enumerator are
+  \begin{description}
+    \item [kUnique] This chain does not share daughter chains with any other through chains
+    \item [kMultipleUS] This chain shares an upstream chain with at least one other through chain
+    \item [kMultipleDS] This chain shares a downstream chain with at least one other through chain
+    \item [kMultipleBoth] This chain shares both its upstream and downstream chains with other through chains
+  \end{description}
+  Note that constellations are possible where e.g. there are 3 chains, one of type \texttt{kMultipleUS}, one of type \texttt{kMultipleDS},
+  and the last of type \texttt{kMultipleBoth}.
+  \item [get\_mapper\_name()] The name of the mapper that last modified this primary chain object.
+  \item [GetUSDaughter()] For a through chain returns the upstream daughter chain.
+  \item [GetDSDaughter()] For a through chain returns the downstream daughter chain.
+  \item [GetMatchedTracks()] Returns a vector of all matched tracks for this chain. There is not a single matched track since as described in section ref{sec:tmprocess},
+        tracks for all feasible PID hypotheses have to be created.
+  \item [GetPIDTrack()] Returns the PID'd track (if one exists in the chain).
+  \item [GetFittedTrack()] Returns the fitted track (if one exists in the chain).
+\end{description}
+
+Note that the global track object also holds references to constituent tracks, so there are two different
+ways one could obtain upstream and downstream tracks from a through primary chain. For example to obtain PID'd tracks:
+\begin{verbatim}
+through_chain->GetPIDTrack()->GetConstituentTracks()
+\end{verbatim}
+or
+\begin{verbatim}
+through_chain->GetUSDaughter()->GetPIDTrack()
+through_chain->GetDSDaughter()->GetPIDTrack()
+\end{verbatim}
+
+\subsection{Identifying Decay Candidates}\label{subsec:decay_candidates}
+
+Primary chains also allow identifying events that may contain a $\pi \rightarrow \mu$ or $\mu \rightarrow e$ decay. If an upstream
+and downstream track come out of Particle ID with different PIDs, no PID'd through track can be created, so an existing through
+chain will return NULL from \texttt{through\_chain->GetPIDTrack()}. In this case, a decay candidate can be identified e.g. by testing that
+\begin{verbatim}
+through_chain->GetUSDaughter()->GetPIDTrack()->get_pid()
+    == MAUS::DataStructure::Global::kPiPlus
+\end{verbatim}
+and
+\begin{verbatim}
+through_chain->GetDSDaughter()->GetPIDTrack()->get_pid()
+    == MAUS::DataStructure::Global::kMuPlus
+\end{verbatim}
+
+Additional candidates could be identified from orphan chains, but caution is advised, as there is a good chance these might be from different particles rather than decays.
+
+\section{Tracks and Space Points from Local Reconstruction}\label{sec:lr_tracks_spacepoints}
+
+While tracks and space points from the local reconstruction which have ended up in matched tracks can be accessed as constituent tracks of matched tracks and
+as space points linked to by the track points in matched tracks, there might be a need to access all tracks and spacepoints from local reconstruction,
+including those that were not successfully matched. Since they can't be associated with a primary chain, they can be accessed direcly from the global event via
+\begin{verbatim}
+std::vector<MAUS::DataStructure::Global::Track*> GetLRTracks()
+std::vector<MAUS::DataStructure::Global::Track*>
+    GetLRTracks(MAUS::DataStructure::Global::DetectorPoint detector)
+\end{verbatim}
+and
+\begin{verbatim}
+std::vector<MAUS::DataStructure::Global::SpacePoint*> GetLRSpacePoints()
+std::vector<MAUS::DataStructure::Global::SpacePoint*>
+    GetLRSpacePoints(MAUS::DataStructure::Global::DetectorPoint detector)
+\end{verbatim}
+where the optional argument detector allows specifying from which detectors tracks or space points should be retrieved.

=== modified file 'doc/maus_user_guide.pdf'
Binary files doc/maus_user_guide.pdf	2017-05-15 14:17:05 +0000 and doc/maus_user_guide.pdf	2017-06-08 17:10:31 +0000 differ
=== modified file 'src/common_cpp/Recon/Global/TrackMatching.cc'
--- src/common_cpp/Recon/Global/TrackMatching.cc	2017-03-17 12:44:22 +0000
+++ src/common_cpp/Recon/Global/TrackMatching.cc	2017-06-08 17:10:31 +0000
@@ -503,7 +503,7 @@
       // If there are multiple matches, this checks if they could have been from the same particle
       AddIfConsistent(temp_spacepoints, hypothesis_track);
     } catch (Exceptions::Exception exc) {
-      Squeak::mout(Squeak::error) << exc.what() << std::endl;
+      Squeak::mout(Squeak::debug) << exc.what() << std::endl;
     }
   }
 }
@@ -528,7 +528,7 @@
       GlobalTools::propagate(x_in, tof1_z - 25.0, field, _max_step_size, pid,
                              _energy_loss);
     } catch (Exceptions::Exception exc) {
-      Squeak::mout(Squeak::error) << exc.what() << std::endl;
+      Squeak::mout(Squeak::debug) << exc.what() << std::endl;
     }
     // Calculate the distance to TOF0 and the approximate distance the particle travelled
     double z_distance = tof1_z - spacepoints.at(0)->get_position().Z();
@@ -626,7 +626,7 @@
             << x_in[2] - first_hit_pos.Y() << std::endl;
       }
     } catch (Exceptions::Exception exc) {
-      Squeak::mout(Squeak::error) << exc.what() << std::endl;
+      Squeak::mout(Squeak::debug) << exc.what() << std::endl;
     }
     if (matched) {
       break;

=== modified file 'src/common_cpp/Recon/SciFi/PatternRecognition.cc'
--- src/common_cpp/Recon/SciFi/PatternRecognition.cc	2017-05-12 14:13:30 +0000
+++ src/common_cpp/Recon/SciFi/PatternRecognition.cc	2017-06-08 17:10:31 +0000
@@ -37,6 +37,8 @@
 
 // MAUS headers
 #include "src/common_cpp/Recon/SciFi/PatternRecognition.hh"
+#include "src/common_cpp/Recon/SciFi/LeastSquaresFitter.hh"
+#include "src/common_cpp/Recon/SciFi/RootFitter.hh"
 #include "src/common_cpp/DataStructure/ThreeVector.hh"
 #include "src/common_cpp/Utils/Globals.hh"
 #include "src/common_cpp/Globals/GlobalsManager.hh"
@@ -67,6 +69,8 @@
                                           _down_helical_pr_on(true),
                                           _sp_search_on(false),
                                           _s_error_method(0),
+                                          _line_fitter(0),
+                                          _circle_fitter(0),
                                           _verb(0),
                                           _n_trackers(2),
                                           _n_stations(5),
@@ -105,6 +109,8 @@
   _down_helical_pr_on = true;
   _sp_search_on = false;
   _s_error_method = 0;
+  _line_fitter = 0;
+  _circle_fitter = 0;
   _verb = 0;
   _n_trackers = 2;
   _n_stations = 5;
@@ -165,6 +171,8 @@
     Json::Value *json = Globals::GetConfigurationCards();
     _sp_search_on = (*json)["SciFiPatRecMissingSpSearchOn"].asBool();
     _s_error_method = (*json)["SciFiPatRecSErrorMethod"].asInt();
+    _line_fitter = (*json)["SciFiPatRecLineFitter"].asInt();
+    _circle_fitter = (*json)["SciFiPatRecCircleFitter"].asInt();
     _verb = (*json)["SciFiPatRecVerbosity"].asInt();
     _n_trackers = (*json)["SciFinTrackers"].asInt();
     _n_stations = (*json)["SciFinStations"].asInt();
@@ -628,8 +636,13 @@
   // Fit track
   SimpleLine line_x, line_y;
   TMatrixD cov_x(2, 2), cov_y(2, 2);
-  LeastSquaresFitter::linear_fit(z, x, x_err, line_x, cov_x);
-  LeastSquaresFitter::linear_fit(z, y, y_err, line_y, cov_y);
+  if (_line_fitter == 0) { // Custom LSQ fitter
+    LeastSquaresFitter::linear_fit(z, x, x_err, line_x, cov_x);
+    LeastSquaresFitter::linear_fit(z, y, y_err, line_y, cov_y);
+  } else if (_line_fitter == 1) { // ROOT based linear fitter
+    RootFitter::FitLineLinear(z, x, x_err, line_x, cov_x);
+    RootFitter::FitLineLinear(z, y, y_err, line_y, cov_y);
+  }
 
   // Squash the covariances of each fit into one vector
   std::vector<double> covariance;
@@ -754,8 +767,24 @@
   TMatrixD cov_circle(3, 3); // The covariance matrix for the circle parameters alpha, beta, gamma
   double s1to4_error = _sd_1to4 * _circle_error_w;
   double s5_error = _sd_5 * _circle_error_w;
-  bool good_radius = LeastSquaresFitter::circle_fit(s1to4_error, s5_error, _R_res_cut,
-                                                    spnts, c_trial, cov_circle);
+
+  bool good_radius = false;
+  if (_circle_fitter == 0) {
+    // With a custom linear least squares fit
+    good_radius = LeastSquaresFitter::circle_fit(s1to4_error, s5_error, _R_res_cut,
+                                                 spnts, c_trial, cov_circle);
+  } else if (_circle_fitter == 1) {
+    // With a ROOT MINUIT based fitter
+    std::vector<double> x;
+    std::vector<double> y;
+    for (auto sp : spnts) {
+        x.push_back(sp->get_position().x());
+        y.push_back(sp->get_position().y());
+    }
+    good_radius = RootFitter::FitCircleMinuit(x, y, c_trial, cov_circle);
+    if (c_trial.get_R() > _R_res_cut)
+        good_radius = false;
+  }
 
   // If the radius calculated is too large or chisq fails, return NULL
   if ( !good_radius || !( c_trial.get_chisq() / ( (2*n_points) - 3 ) < _circle_chisq_cut ) ) {
@@ -884,7 +913,11 @@
   }
 
   // Fit ds and dz to a straight line, to get the gradient, which equals ds/dz
-  LeastSquaresFitter::linear_fit(z_i, s_i, sigma_s, line_sz, cov_sz);
+  if (_line_fitter == 0) { // Custom LSQ fitter
+    LeastSquaresFitter::linear_fit(z_i, s_i, sigma_s, line_sz, cov_sz);
+  } else if (_line_fitter == 1) { // ROOT based linear fitter
+    RootFitter::FitLineLinear(z_i, s_i, sigma_s, line_sz, cov_sz);
+  }
 
   // Check linear fit passes chisq test
   if ( !(line_sz.get_chisq() / ( n_points - 2 ) < _sz_chisq_cut ) ) {

=== modified file 'src/common_cpp/Recon/SciFi/PatternRecognition.hh'
--- src/common_cpp/Recon/SciFi/PatternRecognition.hh	2017-05-12 14:13:30 +0000
+++ src/common_cpp/Recon/SciFi/PatternRecognition.hh	2017-06-08 17:10:31 +0000
@@ -40,7 +40,6 @@
 #include "gtest/gtest_prod.h"
 
 // MAUS headers
-#include "src/common_cpp/Recon/SciFi/LeastSquaresFitter.hh"
 #include "src/common_cpp/Recon/SciFi/SciFiTools.hh"
 #include "src/common_cpp/Recon/SciFi/SimpleLine.hh"
 #include "src/common_cpp/Recon/SciFi/SimpleCircle.hh"
@@ -375,6 +374,18 @@
     /** @brief Set the boolean controlling if we search for missed helical seed spacepoints */
     void set_sp_search_on(bool sp_search_on) { _sp_search_on = sp_search_on; }
 
+    /** @brief Return the line fit method */
+    int get_line_fitter() { return _line_fitter; }
+
+    /** @brief Set line fit method */
+    void set_line_fitter(int line_fitter) { _line_fitter = line_fitter; }
+
+    /** @brief Get the circle fit method */
+    int get_circle_fitter() { return _circle_fitter; }
+
+    /** @brief Set circle fit method */
+    void set_circle_fitter(int circle_fitter) { _circle_fitter = circle_fitter; }
+
     /** @brief Return the verbosity level */
     bool get_verbosity() { return _verb; }
 
@@ -420,6 +431,8 @@
     bool _down_helical_pr_on;   /** Downstream Helical pattern recogntion on or off */
     bool _sp_search_on;         /** Do we seach for seed spoints missed by helical fit? */
     int _s_error_method;        /** How to calc error on s, 0 = station res, 1 = error prop */
+    int _line_fitter;           /** Line fitter, 0 = custom lsq, 1 = ROOT */
+    int _circle_fitter;         /** Circle fitter, 0 = custom lsq, 1 = MINUIT */
     int _verb;                  /** Verbosity: 0=little, 1=more couts */
     int _n_trackers;            /** Number of trackers */
     int _n_stations;            /** Number of stations per tracker */
@@ -440,7 +453,7 @@
     double _Pt_max;             /** MeV/c max Pt for h tracks (given by R_max = 150mm) */
     double _Pz_min;             /** MeV/c min Pz for helical tracks (this is a guess) */
     double _missing_sp_cut;     /** Dist (mm) below which a missing spoint should be added to trk*/
-    // LeastSquaresFitter _lsq;  /** The linear least squares fitting class instance */
+
     TFile* _rfile;   /** A ROOT file pointer for dumping residuals to in debug mode */
     TH1D* _hx;       /** histo of x residuals taken during straight road cut stage */
     TH1D* _hy;       /** histo of y residuals taken during straight road cut stage */

=== added file 'src/common_cpp/Recon/SciFi/RootFitter.cc'
--- src/common_cpp/Recon/SciFi/RootFitter.cc	1970-01-01 00:00:00 +0000
+++ src/common_cpp/Recon/SciFi/RootFitter.cc	2017-06-08 17:10:31 +0000
@@ -0,0 +1,104 @@
+/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
+ *
+ * MAUS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * MAUS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+// ROOT headers
+#include "TLinearFitter.h"
+#include "TMatrixTBase.h"
+#include "TMath.h"
+#include "Math/Functor.h"
+#include "Fit/Fitter.h"
+
+// MAUS headers
+#include "src/common_cpp/Recon/SciFi/RootFitter.hh"
+
+namespace RootFitter {
+
+bool FitLineLinear(const std::vector<double>& x, const std::vector<double>& y,
+                   const std::vector<double>& yerr, MAUS::SimpleLine& line, TMatrixD& cov_matrix) {
+
+  TLinearFitter lf(1);
+  lf.SetFormula("pol1");
+  // Have to copy data into new elements as AssignData wants non-const arguments for some reason
+  std::vector<double> x1 = x;
+  std::vector<double> y1 = y;
+  std::vector<double> yerr1 = yerr;
+  lf.AssignData(x.size(), 1, &x1[0], &y1[0], &yerr1[0]);
+  lf.Eval();
+
+  TVectorD params;
+  TVectorD errors;
+  lf.GetParameters(params);
+  lf.GetErrors(errors);
+
+  line.set_c(params[0]);
+  line.set_m(params[1]);
+  line.set_c_err(errors[0]);
+  line.set_m_err(errors[1]);
+  line.set_chisq(lf.GetChisquare());
+  lf.GetCovarianceMatrix(cov_matrix);
+
+  return true;
+}
+
+bool FitCircleMinuit(const std::vector<double>& x, const std::vector<double>& y,
+                     MAUS::SimpleCircle& circ, TMatrixD& cov_matrix) {
+
+  auto Chi2Function = [&x, &y](const double *par) {
+    // Minimisation function computing the sum of squares of residuals
+    // looping over the points
+    double f = 0.0;
+    for (size_t i = 0; i < x.size(); i++) {
+        double u = x[i] - par[0];
+        double v = y[i] - par[1];
+        double dr = par[2] - std::sqrt(u*u+v*v);
+        f += dr*dr;
+    }
+    return f;
+  };
+
+  // Wrap chi2 function in a function object for the fit
+  // 3 is the number of fit parameters (size of array par)
+  ROOT::Math::Functor fcn(Chi2Function, 3);
+  ROOT::Fit::Fitter fitter;
+  double pStart[3] = {0, 0, 1};
+  fitter.SetFCN(fcn, pStart);
+  fitter.Config().ParSettings(0).SetName("x0");
+  fitter.Config().ParSettings(1).SetName("y0");
+  fitter.Config().ParSettings(2).SetName("R");
+
+  // do the fit
+  bool ok = fitter.FitFCN();
+  if (!ok) {
+    return ok;
+  }
+  const ROOT::Fit::FitResult & result = fitter.Result();
+  result.Print(std::cout);
+
+  // Create a circle object with the results
+  circ = MAUS::SimpleCircle(result.Parameter(0), result.Parameter(1),
+                            result.Parameter(2));
+  circ.set_x0_err(result.Error(0));
+  circ.set_y0_err(result.Error(1));
+  circ.set_R_err(result.Error(2));
+  // circ.set_chisq(result.Chi2());
+  circ.set_chisq(result.MinFcnValue());
+  circ.set_pvalue(result.Prob());
+  result.GetCovarianceMatrix(cov_matrix);
+
+  return true;
+}
+}

=== added file 'src/common_cpp/Recon/SciFi/RootFitter.hh'
--- src/common_cpp/Recon/SciFi/RootFitter.hh	1970-01-01 00:00:00 +0000
+++ src/common_cpp/Recon/SciFi/RootFitter.hh	2017-06-08 17:10:31 +0000
@@ -0,0 +1,57 @@
+/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
+ *
+ * MAUS is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * MAUS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef ROOTFITTER_HH
+#define ROOTFITTER_HH
+
+// C++ headers
+#include <vector>
+
+// ROOT headers
+#include "TMatrixD.h"
+
+// MAUS headers
+#include "src/common_cpp/Recon/SciFi/SimpleLine.hh"
+#include "src/common_cpp/Recon/SciFi/SimpleCircle.hh"
+
+/** @namespace RootFitter
+ *
+ * ROOT fitting routines
+ */
+namespace RootFitter {
+
+/** @brief Fit a straight line using the ROOT TLinearFit class
+  * @param[in] x x coordinates (the indepepndent variable)
+  * @param[in] y y coordinates (the depepndent variable)
+  * @param[in] yerr Errors on the y coordinates
+  * @param[out] line The fit result including errors
+  * @param[out] cov_matrix The 2*2 covariance matrix of the returned parameters
+  */
+bool FitLineLinear(const std::vector<double>& x, const std::vector<double>& y,
+                   const std::vector<double>& yerr, MAUS::SimpleLine& line, TMatrixD& cov_matrix);
+
+/** @brief Fit a circle using the ROOT MINUIT minimiser class
+  * @param[in] x x coordinates (a dependent variable)
+  * @param[in] y y coordinates (a depepndent variable)
+  * @param[out] circ The fit result
+  * @param[out] cov_matrix The 3*3 covariance matrix of the returned parameters
+  */
+bool FitCircleMinuit(const std::vector<double>& x, const std::vector<double>& y,
+                     MAUS::SimpleCircle& circ, TMatrixD& cov_matrix);
+}
+
+#endif

=== modified file 'src/common_cpp/Recon/SciFi/SimpleCircle.cc'
--- src/common_cpp/Recon/SciFi/SimpleCircle.cc	2013-07-09 17:12:03 +0000
+++ src/common_cpp/Recon/SciFi/SimpleCircle.cc	2017-06-08 17:10:31 +0000
@@ -37,6 +37,7 @@
   _kappa_err = 0.0;
   _delta_R = 0.0;
   _chisq = 0.0;
+  _pvalue = 0.0;
 }
 
 SimpleCircle::SimpleCircle(double x0, double y0, double R) {
@@ -56,12 +57,13 @@
   _kappa_err = 0.0;
   _delta_R = 0.0;
   _chisq = 0.0;
+  _pvalue = 0.0;
 }
 
 SimpleCircle::SimpleCircle(double x0, double x0_err, double y0, double y0_err, double R,
                            double R_err, double alpha, double alpha_err, double beta,
                            double beta_err, double gamma, double gamma_err, double kappa,
-                           double kappa_err, double delta_R, double chisq) {
+                           double kappa_err, double delta_R, double chisq, double pvalue) {
   _x0 = x0;
   _x0_err = x0_err;
   _y0 = y0;
@@ -78,6 +80,7 @@
   _kappa_err = kappa_err;
   _delta_R = delta_R;
   _chisq = chisq;
+  _pvalue = pvalue;
 }
 
 // Destructor
@@ -100,12 +103,14 @@
   _kappa_err = 0.0;
   _delta_R = 0.0;
   _chisq = 0.0;
+  _pvalue = 0.0;
 }
 
 void SimpleCircle::set_parameters(double x0, double x0_err, double y0, double y0_err,
                                   double R, double R_err, double alpha, double alpha_err,
                                   double beta, double beta_err, double gamma, double gamma_err,
-                                  double kappa, double kappa_err, double delta_R, double chisq) {
+                                  double kappa, double kappa_err, double delta_R,
+                                  double chisq, double pvalue) {
   _x0 = x0;
   _x0_err = x0_err;
   _y0 = y0;
@@ -122,6 +127,7 @@
   _kappa_err = kappa_err;
   _delta_R = delta_R;
   _chisq = chisq;
+  _pvalue = pvalue;
 }
 
 } // ~namespace MAUS

=== modified file 'src/common_cpp/Recon/SciFi/SimpleCircle.hh'
--- src/common_cpp/Recon/SciFi/SimpleCircle.hh	2013-07-09 17:12:03 +0000
+++ src/common_cpp/Recon/SciFi/SimpleCircle.hh	2017-06-08 17:10:31 +0000
@@ -33,9 +33,10 @@
     SimpleCircle(double x0, double y0, double R);
 
     /** Second parameter constructor */
-    SimpleCircle(double x0, double x0_err, double y0, double y0_err, double R, double R_err,
-                 double alpha, double alpha_err, double beta, double beta_err, double gamma,
-                 double gamma_err, double kappa, double kappa_err, double delta_R, double chisq);
+    SimpleCircle(double x0, double x0_err, double y0, double y0_err, double R,
+                 double R_err, double alpha, double alpha_err, double beta,
+                 double beta_err, double gamma, double gamma_err, double kappa,
+                 double kappa_err, double delta_R, double chisq, double pvalue);
 
     /** Destructor */
     ~SimpleCircle();  // Default destructor
@@ -57,6 +58,7 @@
     double get_kappa_err() const { return _kappa_err; }
 
     double get_chisq() const { return _chisq; }
+    double get_pvalue() const { return _pvalue; }
     double get_delta_R() const { return _delta_R; }
 
     // Setters
@@ -76,12 +78,14 @@
     void set_kappa(double kappa) { _kappa = kappa; }
     void set_kappa_err(double kappa_err) { _kappa_err = kappa_err; }
     void set_chisq(double chisq) { _chisq = chisq; }
+    void set_pvalue(double pvalue) { _pvalue = pvalue; }
     void set_delta_R(double delta_R) { _delta_R = delta_R; }
 
     void set_parameters(double x0, double x0_err, double y0, double y0_err,
                         double R, double R_err, double alpha, double alpha_err,
                         double beta, double beta_err, double gamma, double gamma_err,
-                        double kappa, double kappa_err, double delta_R, double chisq);
+                        double kappa, double kappa_err, double delta_R,
+                        double chisq, double pvalue);
 
   private:
     double _x0;
@@ -100,6 +104,7 @@
     double _kappa_err;
     double _delta_R;
     double _chisq;
+    double _pvalue;
 };
 
 } // ~namespace MAUS

=== modified file 'src/common_py/ConfigurationDefaults.py'
--- src/common_py/ConfigurationDefaults.py	2017-05-11 12:09:56 +0000
+++ src/common_py/ConfigurationDefaults.py	2017-06-08 17:10:31 +0000
@@ -303,6 +303,8 @@
 SciFiPatRecMissingSpCut = 2 # Distance (mm) below which a missing spoint should added to a track
 SciFiPatRecSErrorMethod = 0 # How to calc error on s, 0 = station res, 1 = error prop
 SciFiPatRecVerbosity = 0 # The verbosity of the pat rec (0 - quiet, 1 - more)
+SciFiPatRecLineFitter = 0 # Choose the patrec straight line fitter, 0 = custom lsq, 1 = ROOT
+SciFiPatRecCircleFitter = 0 # Choose the patrec circle fitter, 0 = custom lsq, 1 = MINUIT
 SciFiStraightRoadCut = 7.0 # The road cut in pat rec for straights (mm)
 SciFiStraightChi2Cut = 50.0 # Chi^2 on pat rec straight track fit
 SciFiRadiusResCut = 150.0 # Helix radius cut (mm) for pattern recognition

=== modified file 'tests/cpp_unit/Recon/SciFi/PatternRecognitionTest.cc'
--- tests/cpp_unit/Recon/SciFi/PatternRecognitionTest.cc	2017-05-02 16:26:10 +0000
+++ tests/cpp_unit/Recon/SciFi/PatternRecognitionTest.cc	2017-06-08 17:10:31 +0000
@@ -36,7 +36,7 @@
   virtual void SetUp()    {}
   virtual void TearDown() {}
 
-  std::vector<SciFiSpacePoint*> set_up_spacepoints() {
+  std::vector<SciFiSpacePoint*> set_up_single_track_spacepoints() {
     SciFiSpacePoint *sp1 = new SciFiSpacePoint();
     SciFiSpacePoint *sp2 = new SciFiSpacePoint();
     SciFiSpacePoint *sp3 = new SciFiSpacePoint();
@@ -87,6 +87,139 @@
 
     return spnts;
   }
+
+  std::vector<SciFiSpacePoint*> set_up_multiple_track_spacepoints() {
+    // Here we use the npe just to keep track of which tracker and
+    // track number each sp belongs to
+    std::vector<SciFiSpacePoint*> spnts_t1_trk1;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t1_trk1.push_back(new SciFiSpacePoint());
+      spnts_t1_trk1[i]->set_tracker(0);
+      spnts_t1_trk1[i]->set_station(5-i);
+      spnts_t1_trk1[i]->set_used(false);
+      spnts_t1_trk1[i]->set_npe(11);
+    }
+    std::vector<SciFiSpacePoint*> spnts_t1_trk2;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t1_trk2.push_back(new SciFiSpacePoint());
+      spnts_t1_trk2[i]->set_tracker(0);
+      spnts_t1_trk2[i]->set_station(5-i);
+      spnts_t1_trk2[i]->set_used(false);
+      spnts_t1_trk2[i]->set_npe(12);
+    }
+    std::vector<SciFiSpacePoint*> spnts_t1_trk3;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t1_trk3.push_back(new SciFiSpacePoint());
+      spnts_t1_trk3[i]->set_tracker(0);
+      spnts_t1_trk3[i]->set_station(5-i);
+      spnts_t1_trk3[i]->set_used(false);
+      spnts_t1_trk3[i]->set_npe(13);
+    }
+    std::vector<SciFiSpacePoint*> spnts_t1_trk4;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t1_trk4.push_back(new SciFiSpacePoint());
+      spnts_t1_trk4[i]->set_tracker(0);
+      spnts_t1_trk4[i]->set_station(5-i);
+      spnts_t1_trk4[i]->set_used(false);
+      spnts_t1_trk4[i]->set_npe(14);
+    }
+
+    std::vector<SciFiSpacePoint*> spnts_t2_trk1;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t2_trk1.push_back(new SciFiSpacePoint());
+      spnts_t2_trk1[i]->set_tracker(1);
+      spnts_t2_trk1[i]->set_station(i+1);
+      spnts_t2_trk1[i]->set_used(false);
+      spnts_t2_trk1[i]->set_npe(21);
+    }
+    std::vector<SciFiSpacePoint*> spnts_t2_trk2;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t2_trk2.push_back(new SciFiSpacePoint());
+      spnts_t2_trk2[i]->set_tracker(1);
+      spnts_t2_trk2[i]->set_station(i+1);
+      spnts_t2_trk2[i]->set_used(false);
+      spnts_t2_trk2[i]->set_npe(22);
+    }
+    std::vector<SciFiSpacePoint*> spnts_t2_trk3;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t2_trk3.push_back(new SciFiSpacePoint());
+      spnts_t2_trk3[i]->set_tracker(1);
+      spnts_t2_trk3[i]->set_station(i+1);
+      spnts_t2_trk3[i]->set_used(false);
+      spnts_t2_trk3[i]->set_npe(23);
+    }
+    std::vector<SciFiSpacePoint*> spnts_t2_trk4;
+    for ( size_t i = 0; i < 5; ++i ) {
+      spnts_t2_trk4.push_back(new SciFiSpacePoint());
+      spnts_t2_trk4[i]->set_tracker(1);
+      spnts_t2_trk4[i]->set_station(i+1);
+      spnts_t2_trk4[i]->set_used(false);
+      spnts_t2_trk4[i]->set_npe(24);
+    }
+
+    // Spill 4, mu plus
+    spnts_t1_trk1[0]->set_position(ThreeVector(0.0, 66.44, 1100.0));
+    spnts_t1_trk1[1]->set_position(ThreeVector(-26.4, 47.46, 750.5));
+    spnts_t1_trk1[2]->set_position(ThreeVector(-2.491, 28.47, 450.5));
+    spnts_t1_trk1[3]->set_position(ThreeVector(14.45, 49.18, 200.6));
+    spnts_t1_trk1[4]->set_position(ThreeVector(1.993, 69.03, 0.6523));
+
+    spnts_t2_trk1[0]->set_position(ThreeVector(-26.4, 56.09, 0.6523));
+    spnts_t2_trk1[1]->set_position(ThreeVector(-25.9, -1.726, 200.7));
+    spnts_t2_trk1[2]->set_position(ThreeVector(43.84, -17.26, 450.7));
+    spnts_t2_trk1[3]->set_position(ThreeVector(57.79, 63.85, 750.7));
+    spnts_t2_trk1[4]->set_position(ThreeVector(-32.38, 47.46, 1101));
+
+    // Spill 5, mu plus
+    spnts_t1_trk2[0]->set_position(ThreeVector(-16.44, 15.53, 1100.0));
+    spnts_t1_trk2[1]->set_position(ThreeVector(-19.93, 10.35, 750.5 ));
+    spnts_t1_trk2[2]->set_position(ThreeVector(-15.44, 9.491, 450.5));
+    spnts_t1_trk2[3]->set_position(ThreeVector(-15.44, 12.94, 200.6));
+    spnts_t1_trk2[4]->set_position(ThreeVector(-18.93, 13.81, 0.6523));
+
+    spnts_t2_trk2[0]->set_position(ThreeVector(-4.982, 15.53, 0.6523));
+    spnts_t2_trk2[1]->set_position(ThreeVector(-12.7, -5.609, 200.7));
+    spnts_t2_trk2[2]->set_position(ThreeVector(10.71, -20.28, 450.7));
+    spnts_t2_trk2[3]->set_position(ThreeVector(23.41, 9.491, 750.7));
+    spnts_t2_trk2[4]->set_position(ThreeVector(-12.95, 5.177, 1101 ));
+
+    // Spill 6, mu plus
+    spnts_t1_trk3[0]->set_position(ThreeVector(-50.81, -23.3,  1100));
+    spnts_t1_trk3[1]->set_position(ThreeVector(33.88,  8.628,  750.5));
+    spnts_t1_trk3[2]->set_position(ThreeVector(-41.35, 44.01,  450.5));
+    spnts_t1_trk3[3]->set_position(ThreeVector(-43.84, -31.06, 200.6));
+    spnts_t1_trk3[4]->set_position(ThreeVector(18.93,  -27.61, 0.6523 ));
+
+    spnts_t2_trk3[0]->set_position(ThreeVector(-3.487, 47.46,  0.6523 ));
+    spnts_t2_trk3[1]->set_position(ThreeVector(13.95,  24.16,  200.7));
+    spnts_t2_trk3[2]->set_position(ThreeVector(40.85,  44.87,  450.7));
+    spnts_t2_trk3[3]->set_position(ThreeVector(8.469,  61.26,  750.7 ));
+    spnts_t2_trk3[4]->set_position(ThreeVector(18.43,  26.75,  1101 ));
+
+    // Spill 2, mu minus
+    spnts_t1_trk4[0]->set_position(ThreeVector(-0.4982, 31.06, 1100));
+    spnts_t1_trk4[1]->set_position(ThreeVector(-9.465, -0.8628, 750.5));
+    spnts_t1_trk4[2]->set_position(ThreeVector(20.42, 4.314, 450.5));
+    spnts_t1_trk4[3]->set_position(ThreeVector(11.46, 30.2, 200.6));
+    spnts_t1_trk4[4]->set_position(ThreeVector(-9.465, 25.02, 0.6523));
+
+    spnts_t2_trk4[0]->set_position(ThreeVector(-2.491, -19.85, 0.6523));
+    spnts_t2_trk4[1]->set_position(ThreeVector(12.95, -24.16, 200.7));
+    spnts_t2_trk4[2]->set_position(ThreeVector(18.93, -6.903, 450.7));
+    spnts_t2_trk4[3]->set_position(ThreeVector(-2.491, -4.314, 750.7));
+    spnts_t2_trk4[4]->set_position(ThreeVector(7.971, -25.89, 1101.0));
+
+    std::vector<SciFiSpacePoint*> spnts(spnts_t1_trk1);
+    spnts.insert(spnts.end(), spnts_t1_trk2.begin(), spnts_t1_trk2.end());
+    spnts.insert(spnts.end(), spnts_t1_trk3.begin(), spnts_t1_trk3.end());
+    spnts.insert(spnts.end(), spnts_t1_trk4.begin(), spnts_t1_trk4.end());
+    spnts.insert(spnts.end(), spnts_t2_trk1.begin(), spnts_t2_trk1.end());
+    spnts.insert(spnts.end(), spnts_t2_trk2.begin(), spnts_t2_trk2.end());
+    spnts.insert(spnts.end(), spnts_t2_trk3.begin(), spnts_t2_trk3.end());
+    spnts.insert(spnts.end(), spnts_t2_trk4.begin(), spnts_t2_trk4.end());
+
+    return spnts;
+  }
 };
 
 TEST_F(PatternRecognitionTest, test_constructor) {
@@ -97,6 +230,8 @@
   EXPECT_TRUE(pr._up_helical_pr_on);
   EXPECT_TRUE(pr._down_helical_pr_on);
   EXPECT_FALSE(pr._sp_search_on);
+  EXPECT_EQ(0, pr._line_fitter);
+  EXPECT_EQ(0, pr._circle_fitter);
   EXPECT_EQ(0, pr._verb);
   EXPECT_EQ(2, pr._n_trackers);
   EXPECT_EQ(5, pr._n_stations);
@@ -124,6 +259,8 @@
   EXPECT_TRUE(pr._up_helical_pr_on);
   EXPECT_TRUE(pr._down_helical_pr_on);
   EXPECT_FALSE(pr._sp_search_on);
+  EXPECT_EQ(0, pr._line_fitter);
+  EXPECT_EQ(0, pr._circle_fitter);
   EXPECT_EQ(0, pr._verb);
   EXPECT_EQ(2, pr._n_trackers);
   EXPECT_EQ(5, pr._n_stations);
@@ -149,7 +286,7 @@
   pr.set_parameters_to_default();
 
   // Set up the spacepoints vector
-  std::vector<SciFiSpacePoint*> spnts = set_up_spacepoints();
+  std::vector<SciFiSpacePoint*> spnts = set_up_single_track_spacepoints();
 
   // For a straight fit
   // ------------------
@@ -216,131 +353,14 @@
 }
 */
 
-TEST_F(PatternRecognitionTest, test_multiple_evts_per_trigger) {
+TEST_F(PatternRecognitionTest, test_multiple_evts_per_trigger_lsq) {
 
   PatternRecognition pr;
   pr.set_parameters_to_default();
+  pr.set_circle_fitter(0);
 
   // Set up the spacepoints vector
-  std::vector<SciFiSpacePoint*> spnts_t1_trk1;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t1_trk1.push_back(new SciFiSpacePoint());
-    spnts_t1_trk1[i]->set_tracker(0);
-    spnts_t1_trk1[i]->set_station(5-i);
-    spnts_t1_trk1[i]->set_used(false);
-  }
-  std::vector<SciFiSpacePoint*> spnts_t1_trk2;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t1_trk2.push_back(new SciFiSpacePoint());
-    spnts_t1_trk2[i]->set_tracker(0);
-    spnts_t1_trk2[i]->set_station(5-i);
-    spnts_t1_trk2[i]->set_used(false);
-  }
-  std::vector<SciFiSpacePoint*> spnts_t1_trk3;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t1_trk3.push_back(new SciFiSpacePoint());
-    spnts_t1_trk3[i]->set_tracker(0);
-    spnts_t1_trk3[i]->set_station(5-i);
-    spnts_t1_trk3[i]->set_used(false);
-  }
-  std::vector<SciFiSpacePoint*> spnts_t1_trk4;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t1_trk4.push_back(new SciFiSpacePoint());
-    spnts_t1_trk4[i]->set_tracker(0);
-    spnts_t1_trk4[i]->set_station(5-i);
-    spnts_t1_trk4[i]->set_used(false);
-  }
-
-  std::vector<SciFiSpacePoint*> spnts_t2_trk1;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t2_trk1.push_back(new SciFiSpacePoint());
-    spnts_t2_trk1[i]->set_tracker(1);
-    spnts_t2_trk1[i]->set_station(i+1);
-    spnts_t2_trk1[i]->set_used(false);
-  }
-  std::vector<SciFiSpacePoint*> spnts_t2_trk2;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t2_trk2.push_back(new SciFiSpacePoint());
-    spnts_t2_trk2[i]->set_tracker(1);
-    spnts_t2_trk2[i]->set_station(i+1);
-    spnts_t2_trk2[i]->set_used(false);
-  }
-  std::vector<SciFiSpacePoint*> spnts_t2_trk3;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t2_trk3.push_back(new SciFiSpacePoint());
-    spnts_t2_trk3[i]->set_tracker(1);
-    spnts_t2_trk3[i]->set_station(i+1);
-    spnts_t2_trk3[i]->set_used(false);
-  }
-  std::vector<SciFiSpacePoint*> spnts_t2_trk4;
-  for ( size_t i = 0; i < 5; ++i ) {
-    spnts_t2_trk4.push_back(new SciFiSpacePoint());
-    spnts_t2_trk4[i]->set_tracker(1);
-    spnts_t2_trk4[i]->set_station(i+1);
-    spnts_t2_trk4[i]->set_used(false);
-  }
-
-
-  // Spill 4, mu plus
-  spnts_t1_trk1[0]->set_position(ThreeVector(0.0, 66.44, 1100.0));
-  spnts_t1_trk1[1]->set_position(ThreeVector(-26.4, 47.46, 750.5));
-  spnts_t1_trk1[2]->set_position(ThreeVector(-2.491, 28.47, 450.5));
-  spnts_t1_trk1[3]->set_position(ThreeVector(14.45, 49.18, 200.6));
-  spnts_t1_trk1[4]->set_position(ThreeVector(1.993, 69.03, 0.6523));
-
-  spnts_t2_trk1[0]->set_position(ThreeVector(-26.4, 56.09, 0.6523));
-  spnts_t2_trk1[1]->set_position(ThreeVector(-25.9, -1.726, 200.7));
-  spnts_t2_trk1[2]->set_position(ThreeVector(43.84, -17.26, 450.7));
-  spnts_t2_trk1[3]->set_position(ThreeVector(57.79, 63.85, 750.7));
-  spnts_t2_trk1[4]->set_position(ThreeVector(-32.38, 47.46, 1101));
-
-  // Spill 5, mu plus
-  spnts_t1_trk2[0]->set_position(ThreeVector(-16.44, 15.53, 1100.0));
-  spnts_t1_trk2[1]->set_position(ThreeVector(-19.93, 10.35, 750.5 ));
-  spnts_t1_trk2[2]->set_position(ThreeVector(-15.44, 9.491, 450.5));
-  spnts_t1_trk2[3]->set_position(ThreeVector(-15.44, 12.94, 200.6));
-  spnts_t1_trk2[4]->set_position(ThreeVector(-18.93, 13.81, 0.6523));
-
-  spnts_t2_trk2[0]->set_position(ThreeVector(-4.982, 15.53, 0.6523));
-  spnts_t2_trk2[1]->set_position(ThreeVector(-12.7, -5.609, 200.7));
-  spnts_t2_trk2[2]->set_position(ThreeVector(10.71, -20.28, 450.7));
-  spnts_t2_trk2[3]->set_position(ThreeVector(23.41, 9.491, 750.7));
-  spnts_t2_trk2[4]->set_position(ThreeVector(-12.95, 5.177, 1101 ));
-
-  // Spill 6, mu plus
-  spnts_t1_trk3[0]->set_position(ThreeVector(-50.81, -23.3,  1100));
-  spnts_t1_trk3[1]->set_position(ThreeVector(33.88,  8.628,  750.5));
-  spnts_t1_trk3[2]->set_position(ThreeVector(-41.35, 44.01,  450.5));
-  spnts_t1_trk3[3]->set_position(ThreeVector(-43.84, -31.06, 200.6));
-  spnts_t1_trk3[4]->set_position(ThreeVector(18.93,  -27.61, 0.6523 ));
-
-  spnts_t2_trk3[0]->set_position(ThreeVector(-3.487, 47.46,  0.6523 ));
-  spnts_t2_trk3[1]->set_position(ThreeVector(13.95,  24.16,  200.7));
-  spnts_t2_trk3[2]->set_position(ThreeVector(40.85,  44.87,  450.7));
-  spnts_t2_trk3[3]->set_position(ThreeVector(8.469,  61.26,  750.7 ));
-  spnts_t2_trk3[4]->set_position(ThreeVector(18.43,  26.75,  1101 ));
-
-  // Spill 2, mu minus
-  spnts_t1_trk4[0]->set_position(ThreeVector(-0.4982, 31.06, 1100));
-  spnts_t1_trk4[1]->set_position(ThreeVector(-9.465, -0.8628, 750.5));
-  spnts_t1_trk4[2]->set_position(ThreeVector(20.42, 4.314, 450.5));
-  spnts_t1_trk4[3]->set_position(ThreeVector(11.46, 30.2, 200.6));
-  spnts_t1_trk4[4]->set_position(ThreeVector(-9.465, 25.02, 0.6523));
-
-  spnts_t2_trk4[0]->set_position(ThreeVector(-2.491, -19.85, 0.6523));
-  spnts_t2_trk4[1]->set_position(ThreeVector(12.95, -24.16, 200.7));
-  spnts_t2_trk4[2]->set_position(ThreeVector(18.93, -6.903, 450.7));
-  spnts_t2_trk4[3]->set_position(ThreeVector(-2.491, -4.314, 750.7));
-  spnts_t2_trk4[4]->set_position(ThreeVector(7.971, -25.89, 1101.0));
-
-  std::vector<SciFiSpacePoint*> spnts(spnts_t1_trk1);
-  spnts.insert(spnts.end(), spnts_t1_trk2.begin(), spnts_t1_trk2.end());
-  spnts.insert(spnts.end(), spnts_t1_trk3.begin(), spnts_t1_trk3.end());
-  spnts.insert(spnts.end(), spnts_t1_trk4.begin(), spnts_t1_trk4.end());
-  spnts.insert(spnts.end(), spnts_t2_trk1.begin(), spnts_t2_trk1.end());
-  spnts.insert(spnts.end(), spnts_t2_trk2.begin(), spnts_t2_trk2.end());
-  spnts.insert(spnts.end(), spnts_t2_trk3.begin(), spnts_t2_trk3.end());
-  spnts.insert(spnts.end(), spnts_t2_trk4.begin(), spnts_t2_trk4.end());
+  std::vector<SciFiSpacePoint*> spnts = set_up_multiple_track_spacepoints();
   SciFiEvent evt1;
   evt1.set_spacepoints(spnts);
 
@@ -377,6 +397,14 @@
   EXPECT_EQ(5, htrks[5]->get_num_points());
   EXPECT_EQ(5, htrks[6]->get_num_points());
   EXPECT_EQ(5, htrks[7]->get_num_points());
+  // Check npe, which we used to encode which sp belongs to which tracks
+  for (SciFiHelicalPRTrack* trk : htrks) {
+    std::vector<SciFiSpacePoint*> spnts = trk->get_spacepoints_pointers();
+    EXPECT_NEAR(spnts[1]->get_npe(), spnts[0]->get_npe(), 0.01);
+    EXPECT_NEAR(spnts[2]->get_npe(), spnts[0]->get_npe(), 0.01);
+    EXPECT_NEAR(spnts[3]->get_npe(), spnts[0]->get_npe(), 0.01);
+    EXPECT_NEAR(spnts[4]->get_npe(), spnts[0]->get_npe(), 0.01);
+  }
   EXPECT_NEAR(htrks[0]->get_dsdz(), -0.1156, 0.001);
   EXPECT_NEAR(htrks[1]->get_dsdz(), -0.342, 0.01);
   EXPECT_NEAR(htrks[2]->get_dsdz(), -0.01834, 0.01);
@@ -385,14 +413,74 @@
   EXPECT_NEAR(htrks[5]->get_dsdz(), 0.3126, 0.001);
   EXPECT_NEAR(htrks[6]->get_dsdz(), 0.08396, 0.001);
   EXPECT_NEAR(htrks[7]->get_dsdz(), 0.1257, 0.001);
-
-  // evt descoping will delete the spacepoints
+}
+
+TEST_F(PatternRecognitionTest, test_multiple_evts_per_trigger_rootfit) {
+
+  PatternRecognition pr;
+  pr.set_parameters_to_default();
+  pr.set_circle_fitter(1);
+
+  // Set up the spacepoints vector
+  std::vector<SciFiSpacePoint*> spnts = set_up_multiple_track_spacepoints();
+  SciFiEvent evt1;
+  evt1.set_spacepoints(spnts);
+
+  // Randomise things a bit to make it harder
+  SciFiSpacePoint *sp1, *sp2;
+  sp1 = spnts[3];
+  sp2 = spnts[14];
+  spnts[3] = sp2;
+  spnts[14] = sp1;
+  sp1 = spnts[4];
+  sp2 = spnts[17];
+  spnts[4] = sp2;
+  spnts[17] = sp1;
+
+  // Perform the recon
+  pr.set_up_helical_pr_on(true);
+  pr.set_down_helical_pr_on(true);
+  pr.set_up_straight_pr_on(false);
+  pr.set_down_straight_pr_on(false);
+  pr.process(evt1);
+
+  std::vector<SciFiStraightPRTrack*> strks;
+  std::vector<SciFiHelicalPRTrack*> htrks;
+  strks = evt1.straightprtracks();
+  htrks = evt1.helicalprtracks();
+
+  ASSERT_EQ(8u, htrks.size());
+  EXPECT_EQ(0u, strks.size());
+  EXPECT_EQ(5, htrks[0]->get_num_points());
+  EXPECT_EQ(5, htrks[1]->get_num_points());
+  EXPECT_EQ(5, htrks[2]->get_num_points());
+  EXPECT_EQ(5, htrks[3]->get_num_points());
+  EXPECT_EQ(5, htrks[4]->get_num_points());
+  EXPECT_EQ(5, htrks[5]->get_num_points());
+  EXPECT_EQ(5, htrks[6]->get_num_points());
+  EXPECT_EQ(5, htrks[7]->get_num_points());
+  // Check npe, which we used to encode which sp belongs to which tracks
+  for (SciFiHelicalPRTrack* trk : htrks) {
+    std::vector<SciFiSpacePoint*> spnts = trk->get_spacepoints_pointers();
+    EXPECT_NEAR(spnts[1]->get_npe(), spnts[0]->get_npe(), 0.01);
+    EXPECT_NEAR(spnts[2]->get_npe(), spnts[0]->get_npe(), 0.01);
+    EXPECT_NEAR(spnts[3]->get_npe(), spnts[0]->get_npe(), 0.01);
+    EXPECT_NEAR(spnts[4]->get_npe(), spnts[0]->get_npe(), 0.01);
+  }
+  EXPECT_NEAR(htrks[0]->get_dsdz(), -0.342, 0.01);
+  EXPECT_NEAR(htrks[1]->get_dsdz(), -0.1156, 0.005);
+  EXPECT_NEAR(htrks[2]->get_dsdz(), -0.01834, 0.01);
+  EXPECT_NEAR(htrks[3]->get_dsdz(), -0.1178, 0.01);
+  EXPECT_NEAR(htrks[4]->get_dsdz(), 0.08396, 0.001);
+  EXPECT_NEAR(htrks[5]->get_dsdz(), 0.3126, 0.001);
+  EXPECT_NEAR(htrks[6]->get_dsdz(), 0.1257, 0.001);
+  EXPECT_NEAR(htrks[7]->get_dsdz(), 0.1504, 0.001);
 }
 
 TEST_F(PatternRecognitionTest, test_make_tracks) {
 
   // Set up the spacepoints vector
-  std::vector<SciFiSpacePoint*> spnts_all = set_up_spacepoints();
+  std::vector<SciFiSpacePoint*> spnts_all = set_up_single_track_spacepoints();
   std::vector<SciFiSpacePoint*> spnts;
   spnts.push_back(spnts_all[4]);
   spnts.push_back(spnts_all[1]);
@@ -517,7 +605,7 @@
 TEST_F(PatternRecognitionTest, test_make_4pt_tracks) {
 
   // Set up the spacepoints vector
-  std::vector<SciFiSpacePoint*> spnts_all = set_up_spacepoints();
+  std::vector<SciFiSpacePoint*> spnts_all = set_up_single_track_spacepoints();
   std::vector<SciFiSpacePoint*> spnts;
   spnts.push_back(spnts_all[4]);
   spnts.push_back(spnts_all[2]);
@@ -617,7 +705,7 @@
 TEST_F(PatternRecognitionTest, test_make_3pt_tracks) {
 
   // Set up the spacepoints vector
-  std::vector<SciFiSpacePoint*> spnts_all = set_up_spacepoints();
+  std::vector<SciFiSpacePoint*> spnts_all = set_up_single_track_spacepoints();
   std::vector<SciFiSpacePoint*> spnts;
   spnts.push_back(spnts_all[4]);
   spnts.push_back(spnts_all[1]);
@@ -737,15 +825,75 @@
   }
 }
 
-TEST_F(PatternRecognitionTest, test_make_straight_tracks) {
-
-  int n_stations = 5;
-  int tracker_num = 0;
-  PatternRecognition pr;
-  pr.set_parameters_to_default();
-
-  // Set up the spacepoints vector
-  std::vector<SciFiSpacePoint*> spnts = set_up_spacepoints();
+TEST_F(PatternRecognitionTest, test_make_straight_tracks_lsq) {
+
+  int n_stations = 5;
+  int tracker_num = 0;
+  PatternRecognition pr;
+  pr.set_parameters_to_default();
+  pr.set_line_fitter(0);
+
+  // Set up the spacepoints vector
+  std::vector<SciFiSpacePoint*> spnts = set_up_single_track_spacepoints();
+
+  // Set up the spacepoints by station 2D vector
+  std::vector< std::vector<SciFiSpacePoint*> > spnts_by_station(n_stations);
+  SciFiTools::sort_by_station(spnts, spnts_by_station);
+
+  // Check the spacepoints have setup correctly
+  EXPECT_EQ(spnts[0], spnts_by_station[0][0]);
+  EXPECT_EQ(spnts[1], spnts_by_station[1][0]);
+  EXPECT_EQ(spnts[2], spnts_by_station[2][0]);
+  EXPECT_EQ(spnts[3], spnts_by_station[3][0]);
+  EXPECT_EQ(spnts[4], spnts_by_station[4][0]);
+  EXPECT_EQ(-68.24883333333334, spnts_by_station[0][0]->get_position().x());
+
+  std::vector<int> ignore_stations;
+  std::vector<SciFiStraightPRTrack*> strks;
+
+  // The track parameters that should be reconstructed from the spacepoints
+  int num_points = 5;
+  double x_chisq = 22.87148204;
+  double y_chisq = 20.99052559;
+  double y0 = -58.85201389;
+  double x0 = -68.94108927;
+  double my = 0.03755825;
+  double mx = -0.02902014;
+
+  // Make the track from the spacepoints
+  pr.make_straight_tracks(num_points, tracker_num, ignore_stations, spnts_by_station, strks);
+
+  // Check it matches to within a tolerance epsilon
+  double epsilon = 0.000001;
+  EXPECT_EQ(1u, strks.size());
+  EXPECT_NEAR(x0, strks[0]->get_x0(), epsilon);
+  EXPECT_NEAR(mx, strks[0]->get_mx(), epsilon);
+  EXPECT_NEAR(x_chisq, strks[0]->get_x_chisq(), epsilon);
+  EXPECT_NEAR(y0, strks[0]->get_y0(), epsilon);
+  EXPECT_NEAR(my, strks[0]->get_my(), epsilon);
+  EXPECT_NEAR(y_chisq, strks[0]->get_y_chisq(), epsilon);
+
+  // Tidy up
+  std::vector<SciFiSpacePoint*>::iterator it;
+  for (it = spnts.begin(); it != spnts.end(); ++it) {
+    delete (*it);
+  }
+  std::vector<SciFiStraightPRTrack*>::iterator strack;
+  for (strack = strks.begin(); strack != strks.end(); ++strack) {
+    delete (*strack);
+  }
+}
+
+TEST_F(PatternRecognitionTest, test_make_straight_tracks_rootfit) {
+
+  int n_stations = 5;
+  int tracker_num = 0;
+  PatternRecognition pr;
+  pr.set_parameters_to_default();
+  pr.set_line_fitter(1);
+
+  // Set up the spacepoints vector
+  std::vector<SciFiSpacePoint*> spnts = set_up_single_track_spacepoints();
 
   // Set up the spacepoints by station 2D vector
   std::vector< std::vector<SciFiSpacePoint*> > spnts_by_station(n_stations);

=== modified file 'tests/cpp_unit/Recon/SciFi/SimpleCircleTest.cc'
--- tests/cpp_unit/Recon/SciFi/SimpleCircleTest.cc	2013-07-09 17:12:03 +0000
+++ tests/cpp_unit/Recon/SciFi/SimpleCircleTest.cc	2017-06-08 17:10:31 +0000
@@ -92,9 +92,10 @@
   double kappa_err = 0.8;
   double delta_R = 9;
   double chisq = 10;
+  double pvalue = 11.0;
 
   SimpleCircle circ1(x0, x0_err, y0, y0_err, R, R_err, alpha, alpha_err, beta, beta_err,
-                     gamma, gamma_err, kappa, kappa_err, delta_R, chisq);
+                     gamma, gamma_err, kappa, kappa_err, delta_R, chisq, pvalue);
 
   EXPECT_EQ(circ1.get_x0(), x0);
   EXPECT_EQ(circ1.get_x0_err(), x0_err);
@@ -112,6 +113,7 @@
   EXPECT_EQ(circ1.get_kappa_err(), kappa_err);
   EXPECT_EQ(circ1.get_delta_R(), delta_R);
   EXPECT_EQ(circ1.get_chisq(), chisq);
+  EXPECT_EQ(circ1.get_pvalue(), pvalue);
 }
 
 TEST_F(SimpleCircleTestDS, test_getters_setters_clear) {
@@ -131,6 +133,7 @@
   double kappa_err = 0.8;
   double delta_R = 9.0;
   double chisq = 10.0;
+  double pvalue = 11.0;
 
   SimpleCircle circ1;
 
@@ -150,6 +153,7 @@
   circ1.set_kappa_err(kappa_err);
   circ1.set_delta_R(delta_R);
   circ1.set_chisq(chisq);
+  circ1.set_pvalue(pvalue);
 
   EXPECT_EQ(circ1.get_x0(), x0);
   EXPECT_EQ(circ1.get_x0_err(), x0_err);
@@ -167,6 +171,7 @@
   EXPECT_EQ(circ1.get_kappa_err(), kappa_err);
   EXPECT_EQ(circ1.get_delta_R(), delta_R);
   EXPECT_EQ(circ1.get_chisq(), chisq);
+  EXPECT_EQ(circ1.get_pvalue(), pvalue);
 
   circ1.clear();
 
@@ -186,9 +191,10 @@
   EXPECT_EQ(circ1.get_kappa_err(), 0.0);
   EXPECT_EQ(circ1.get_delta_R(), 0.0);
   EXPECT_EQ(circ1.get_chisq(), 0.0);
+  EXPECT_EQ(circ1.get_pvalue(), 0.0);
 
   circ1.set_parameters(x0, x0_err, y0, y0_err, R, R_err, alpha, alpha_err, beta, beta_err,
-                       gamma, gamma_err, kappa, kappa_err, delta_R, chisq);
+                       gamma, gamma_err, kappa, kappa_err, delta_R, chisq, pvalue);
 
   EXPECT_EQ(circ1.get_x0(), x0);
   EXPECT_EQ(circ1.get_x0_err(), x0_err);
@@ -206,6 +212,7 @@
   EXPECT_EQ(circ1.get_kappa_err(), kappa_err);
   EXPECT_EQ(circ1.get_delta_R(), delta_R);
   EXPECT_EQ(circ1.get_chisq(), chisq);
+  EXPECT_EQ(circ1.get_pvalue(), pvalue);
 }
 
 } // ~namespace MAUS

=== modified file 'tests/py_unit/test_maus_cpp/test_optics/test_optics_model.py'
--- tests/py_unit/test_maus_cpp/test_optics/test_optics_model.py	2015-07-08 10:26:19 +0000
+++ tests/py_unit/test_maus_cpp/test_optics/test_optics_model.py	2017-06-08 17:10:31 +0000
@@ -1,3 +1,5 @@
+#!/usr/bin/env python
+
 #  This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
 # 
 #  MAUS is free software: you can redistribute it and/or modify
@@ -19,20 +21,171 @@
 Test maus_cpp.optics_model
 """
 
+import json
+import unittest
 import os
-import subprocess
-import unittest
+
+import Configuration
+import maus_cpp.globals
+import maus_cpp.mice_module
+import maus_cpp.covariance_matrix
+from maus_cpp.covariance_matrix import CovarianceMatrix
+import maus_cpp.phase_space_vector
+from maus_cpp.phase_space_vector import PhaseSpaceVector
+import maus_cpp.optics_model
+from maus_cpp.optics_model import OpticsModel
+
+REF = {
+    "position":{"x":0., "y":0., "z":1000.001},
+    "momentum":{"x":0., "y":0., "z":1.},
+    "particle_id":-13, "energy":226., "random_seed":0, "time":0.
+}
 
 class OpticsModelTestCase(unittest.TestCase): # pylint: disable=R0904
     """Test maus_cpp.optics_model"""
-    def test_optics_model(self): #pylint: disable=R0201
-        """
-        Test optics model runs in a subprocess to prevent pollution of globals
-        """
-        test = os.path.expandvars(\
+    def setUp(self):
+        """import datacards"""
+        self.test_config = ""
+        self.good_mod = os.path.expandvars(\
               "${MAUS_ROOT_DIR}/tests/py_unit/test_maus_cpp/"+\
-              "test_optics/test_optics_model_subprocess")
-        subprocess.check_call(test)
+              "test_optics/optics_model.dat")
+        self.no_planes_mod = os.path.expandvars("Test.dat")
+
+    def _set_geometry(self, geom_filename):
+        """Geometry setup"""
+        if maus_cpp.globals.has_instance():
+            maus_cpp.globals.death()
+        self.test_config = Configuration.Configuration().getConfigJSON()
+        json_config = json.loads(self.test_config)
+        json_config["simulation_geometry_filename"] = geom_filename
+        json_config["simulation_reference_particle"] = REF
+        json_config["physics_processes"] = "none"
+        self.test_config = json.dumps(json_config)
+        maus_cpp.globals.birth(self.test_config)
+        maus_cpp.globals.set_monte_carlo_mice_modules(
+                                maus_cpp.mice_module.MiceModule(geom_filename))
+
+    def test_init_no_globals(self):
+        """Test maus_cpp.optics_model.__init__() and deallocation"""
+        if maus_cpp.globals.has_instance():
+            maus_cpp.globals.death()
+        try:
+            OpticsModel()
+            self.assertTrue(False, msg="Should throw an exception if globals "+\
+                                       "is not initialised")
+        except RuntimeError:
+            pass
+
+    def test_init_all_okay(self):
+        """Test maus_cpp.optics_model.__init__() and deallocation"""
+        self._set_geometry(self.good_mod)
+        optics = OpticsModel()
+        optics.__init__() # legal, should reinitialise
+
+    def test_transport_covariance_matrix_no_planes(self):
+        """
+        Test maus_cpp.optics_model.Optics().transport_covariance_matrix()
+        with no virtuals
+        """
+        self._set_geometry(self.no_planes_mod)
+        optics = OpticsModel()
+        cm_in = CovarianceMatrix()
+        try:
+            optics.transport_covariance_matrix(cm_in, 2000.)
+            self.assertTrue(False, "Should throw when no virtuals")
+        except RuntimeError:
+            pass
+
+    def test_transport_covariance_matrix_bad_type(self):
+        """
+        Test maus_cpp.optics_model.Optics().transport_covariance_matrix()
+        with non cm
+        """
+        self._set_geometry(self.good_mod)
+        optics = OpticsModel()
+        try:
+            optics.transport_covariance_matrix("should be a cm", 2000.)
+            self.assertTrue(False, "Should throw when wrong type passed")
+        except TypeError:
+            pass
+
+    def test_transport_covariance_matrix(self):
+        """Test maus_cpp.optics_model.Optics().transport_covariance_matrix()"""
+        self._set_geometry(self.good_mod)
+        optics = OpticsModel()
+        cm_in = maus_cpp.covariance_matrix.create_from_penn_parameters(
+                  mass=105.658, momentum=200., emittance_t=6., beta_t=333.,
+                 emittance_l=1., beta_l=10., bz=0.)
+        # check energy, px RMS does not change (no fields)
+        cm_out_1 = optics.transport_covariance_matrix(cm_in, 2000.)
+        for i in [2, 4, 6]:
+            self.assertAlmostEqual(cm_in.get_element(i, i),
+                                  cm_out_1.get_element(i, i), 1,
+                                msg="\nIN\n"+str(cm_in)+"\nOUT\n"+str(cm_out_1))
+        cm_out_2 = optics.transport_covariance_matrix(cm_in, 3000.)
+        for i in range(3):
+            self.assertGreater(abs(cm_out_1.get_element(2*i+1, 2*i+2)), 1.)
+        # check that we have some growth of correlation between e.g. x, px
+        self.assertAlmostEqual(cm_out_1.get_element(3, 4),
+                               cm_out_1.get_element(5, 6), 2)
+        for i in [2, 4, 6]:
+            self.assertAlmostEqual(cm_in.get_element(i, i),
+                                  cm_out_2.get_element(i, i), 1,
+                                msg="\nIN\n"+str(cm_in)+"\nOUT\n"+str(cm_out_2))
+        # check no coupling between phase space planes
+        for i in [2, 4, 5]:
+            self.assertAlmostEqual(
+                  cm_in.get_element(i+1, i)-cm_out_1.get_element(i+1, i), 
+                  cm_out_1.get_element(i+1, i)-cm_out_2.get_element(i+1, i),
+                  1,
+                  msg="\nIN\n"+str(cm_in)+"\nOUT\n"+str(cm_out_2))
+
+
+    def test_transport_phase_space_vector_no_planes(self):
+        """
+        Test maus_cpp.optics_model.Optics().transport_phase_space_vector() with
+        no virtuals
+        """
+        self._set_geometry(self.no_planes_mod)
+        optics = OpticsModel()
+        psv_in = PhaseSpaceVector()
+        try:
+            optics.transport_phase_space_vector(psv_in, 2000.)
+            self.assertTrue(False, "Should throw when no virtuals")
+        except RuntimeError:
+            pass
+
+    def test_transport_phase_space_vector_bad_type(self):
+        """
+        Test maus_cpp.optics_model.Optics().transport_phase_space_vector() with
+        no virtuals
+        """
+        self._set_geometry(self.good_mod)
+        optics = OpticsModel()
+        try:
+            optics.transport_phase_space_vector("should be a psv", 2000.)
+            self.assertTrue(False, "Should throw when wrong type passed")
+        except TypeError:
+            pass
+
+    def test_transport_phase_space_vector(self):
+        """Test maus_cpp.optics_model.Optics().transport_phase_space_vector()"""
+        self._set_geometry(self.good_mod)
+        optics = OpticsModel()
+        psv_in = maus_cpp.phase_space_vector.create_from_coordinates \
+                                                    (0.1, 226., 1., 2., 3., 4.)
+        psv_out = optics.transport_phase_space_vector(psv_in, 2000.)
+        pz = (226.**2-105.6583715**2.-2.**2-4.**2)**0.5
+        dz = 1000.
+        t_expected = psv_in.get_t()+psv_out.get_energy()/pz/300.*dz # c=300.
+        x_expected = psv_in.get_x()+psv_in.get_px()/pz*dz
+        y_expected = psv_in.get_y()+psv_in.get_py()/pz*dz
+        self.assertAlmostEqual(psv_out.get_t()/t_expected, 1, 3)
+        self.assertAlmostEqual(psv_out.get_energy()/psv_in.get_energy(), 1., 5)
+        self.assertAlmostEqual(psv_out.get_x()/x_expected, 1., 3)
+        self.assertAlmostEqual(psv_out.get_px()/psv_in.get_px(), 1., 5)
+        self.assertAlmostEqual(psv_out.get_y()/y_expected, 1., 3)
+        self.assertAlmostEqual(psv_out.get_py()/psv_in.get_py(), 1., 5)
 
 if __name__ == "__main__":
     unittest.main()

=== removed file 'tests/py_unit/test_maus_cpp/test_optics/test_optics_model_subprocess'
--- tests/py_unit/test_maus_cpp/test_optics/test_optics_model_subprocess	2015-07-08 10:06:34 +0000
+++ tests/py_unit/test_maus_cpp/test_optics/test_optics_model_subprocess	1970-01-01 00:00:00 +0000
@@ -1,192 +0,0 @@
-#!/usr/bin/env python
-
-#  This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
-# 
-#  MAUS is free software: you can redistribute it and/or modify
-#  it under the terms of the GNU General Public License as published by
-#  the Free Software Foundation, either version 3 of the License, or
-#  (at your option) any later version.
-# 
-#  MAUS is distributed in the hope that it will be useful,
-#  but WITHOUT ANY WARRANTY; without even the implied warranty of
-#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-#  GNU General Public License for more details.
-# 
-#  You should have received a copy of the GNU General Public License
-#  along with MAUS.  If not, see <http://www.gnu.org/licenses/>.
-
-# pylint: disable=C0103
-
-"""
-Test maus_cpp.optics_model
-"""
-
-import json
-import unittest
-import os
-
-import Configuration
-import maus_cpp.globals
-import maus_cpp.mice_module
-import maus_cpp.covariance_matrix
-from maus_cpp.covariance_matrix import CovarianceMatrix
-import maus_cpp.phase_space_vector
-from maus_cpp.phase_space_vector import PhaseSpaceVector
-import maus_cpp.optics_model
-from maus_cpp.optics_model import OpticsModel
-
-REF = {
-    "position":{"x":0., "y":0., "z":1000.001},
-    "momentum":{"x":0., "y":0., "z":1.},
-    "particle_id":-13, "energy":226., "random_seed":0, "time":0.
-}
-
-class OpticsModelTestCase(unittest.TestCase): # pylint: disable=R0904
-    """Test maus_cpp.optics_model"""
-    def setUp(self):
-        """import datacards"""
-        self.test_config = ""
-        self.good_mod = os.path.expandvars(\
-              "${MAUS_ROOT_DIR}/tests/py_unit/test_maus_cpp/"+\
-              "test_optics/optics_model.dat")
-        self.no_planes_mod = os.path.expandvars("Test.dat")
-
-    def _set_geometry(self, geom_filename):
-        """Geometry setup"""
-        if maus_cpp.globals.has_instance():
-            maus_cpp.globals.death()
-        self.test_config = Configuration.Configuration().getConfigJSON()
-        json_config = json.loads(self.test_config)
-        json_config["simulation_geometry_filename"] = geom_filename
-        json_config["simulation_reference_particle"] = REF
-        json_config["physics_processes"] = "none"
-        self.test_config = json.dumps(json_config)
-        maus_cpp.globals.birth(self.test_config)
-        maus_cpp.globals.set_monte_carlo_mice_modules(
-                                maus_cpp.mice_module.MiceModule(geom_filename))
-
-    def test_init_no_globals(self):
-        """Test maus_cpp.optics_model.__init__() and deallocation"""
-        if maus_cpp.globals.has_instance():
-            maus_cpp.globals.death()
-        try:
-            OpticsModel()
-            self.assertTrue(False, msg="Should throw an exception if globals "+\
-                                       "is not initialised")
-        except RuntimeError:
-            pass
-
-    def test_init_all_okay(self):
-        """Test maus_cpp.optics_model.__init__() and deallocation"""
-        self._set_geometry(self.good_mod)
-        optics = OpticsModel()
-        optics.__init__() # legal, should reinitialise
-
-    def test_transport_covariance_matrix_no_planes(self):
-        """
-        Test maus_cpp.optics_model.Optics().transport_covariance_matrix()
-        with no virtuals
-        """
-        self._set_geometry(self.no_planes_mod)
-        optics = OpticsModel()
-        cm_in = CovarianceMatrix()
-        try:
-            optics.transport_covariance_matrix(cm_in, 2000.)
-            self.assertTrue(False, "Should throw when no virtuals")
-        except RuntimeError:
-            pass
-
-    def test_transport_covariance_matrix_bad_type(self):
-        """
-        Test maus_cpp.optics_model.Optics().transport_covariance_matrix()
-        with non cm
-        """
-        self._set_geometry(self.good_mod)
-        optics = OpticsModel()
-        try:
-            optics.transport_covariance_matrix("should be a cm", 2000.)
-            self.assertTrue(False, "Should throw when wrong type passed")
-        except TypeError:
-            pass
-
-    def test_transport_covariance_matrix(self):
-        """Test maus_cpp.optics_model.Optics().transport_covariance_matrix()"""
-        self._set_geometry(self.good_mod)
-        optics = OpticsModel()
-        cm_in = maus_cpp.covariance_matrix.create_from_penn_parameters(
-                  mass=105.658, momentum=200., emittance_t=6., beta_t=333.,
-                 emittance_l=1., beta_l=10., bz=0.)
-        # check energy, px RMS does not change (no fields)
-        cm_out_1 = optics.transport_covariance_matrix(cm_in, 2000.)
-        for i in [2, 4, 6]:
-            self.assertAlmostEqual(cm_in.get_element(i, i),
-                                  cm_out_1.get_element(i, i), 1,
-                                msg="\nIN\n"+str(cm_in)+"\nOUT\n"+str(cm_out_1))
-        cm_out_2 = optics.transport_covariance_matrix(cm_in, 3000.)
-        for i in range(3):
-            self.assertGreater(abs(cm_out_1.get_element(2*i+1, 2*i+2)), 1.)
-        # check that we have some growth of correlation between e.g. x, px
-        self.assertAlmostEqual(cm_out_1.get_element(3, 4),
-                               cm_out_1.get_element(5, 6), 2)
-        for i in [2, 4, 6]:
-            self.assertAlmostEqual(cm_in.get_element(i, i),
-                                  cm_out_2.get_element(i, i), 1,
-                                msg="\nIN\n"+str(cm_in)+"\nOUT\n"+str(cm_out_2))
-        # check no coupling between phase space planes
-        for i in [2, 4, 5]:
-            self.assertAlmostEqual(
-                  cm_in.get_element(i+1, i)-cm_out_1.get_element(i+1, i), 
-                  cm_out_1.get_element(i+1, i)-cm_out_2.get_element(i+1, i),
-                  1,
-                  msg="\nIN\n"+str(cm_in)+"\nOUT\n"+str(cm_out_2))
-
-
-    def test_transport_phase_space_vector_no_planes(self):
-        """
-        Test maus_cpp.optics_model.Optics().transport_phase_space_vector() with
-        no virtuals
-        """
-        self._set_geometry(self.no_planes_mod)
-        optics = OpticsModel()
-        psv_in = PhaseSpaceVector()
-        try:
-            optics.transport_phase_space_vector(psv_in, 2000.)
-            self.assertTrue(False, "Should throw when no virtuals")
-        except RuntimeError:
-            pass
-
-    def test_transport_phase_space_vector_bad_type(self):
-        """
-        Test maus_cpp.optics_model.Optics().transport_phase_space_vector() with
-        no virtuals
-        """
-        self._set_geometry(self.good_mod)
-        optics = OpticsModel()
-        try:
-            optics.transport_phase_space_vector("should be a psv", 2000.)
-            self.assertTrue(False, "Should throw when wrong type passed")
-        except TypeError:
-            pass
-
-    def test_transport_phase_space_vector(self):
-        """Test maus_cpp.optics_model.Optics().transport_phase_space_vector()"""
-        self._set_geometry(self.good_mod)
-        optics = OpticsModel()
-        psv_in = maus_cpp.phase_space_vector.create_from_coordinates \
-                                                    (0.1, 226., 1., 2., 3., 4.)
-        psv_out = optics.transport_phase_space_vector(psv_in, 2000.)
-        pz = (226.**2-105.6583715**2.-2.**2-4.**2)**0.5
-        dz = 1000.
-        t_expected = psv_in.get_t()+psv_out.get_energy()/pz/300.*dz # c=300.
-        x_expected = psv_in.get_x()+psv_in.get_px()/pz*dz
-        y_expected = psv_in.get_y()+psv_in.get_py()/pz*dz
-        self.assertAlmostEqual(psv_out.get_t()/t_expected, 1, 3)
-        self.assertAlmostEqual(psv_out.get_energy()/psv_in.get_energy(), 1., 5)
-        self.assertAlmostEqual(psv_out.get_x()/x_expected, 1., 3)
-        self.assertAlmostEqual(psv_out.get_px()/psv_in.get_px(), 1., 5)
-        self.assertAlmostEqual(psv_out.get_y()/y_expected, 1., 3)
-        self.assertAlmostEqual(psv_out.get_py()/psv_in.get_py(), 1., 5)
-
-if __name__ == "__main__":
-    unittest.main()
-

=== modified file 'tests/unit_tests.bash'
--- tests/unit_tests.bash	2017-02-20 14:14:11 +0000
+++ tests/unit_tests.bash	2017-06-08 17:10:31 +0000
@@ -41,7 +41,14 @@
     maus_coverage_extra="--with-coverage --cover-html --cover-html-dir=${MAUS_ROOT_DIR}/doc/python_coverage --cover-inclusive"
 fi
 
-python ${MAUS_THIRD_PARTY}/third_party/install/bin/nosetests -v ${MAUS_ROOT_DIR}/build $maus_coverage_extra
+python ${MAUS_THIRD_PARTY}/third_party/install/bin/nosetests -v ${MAUS_ROOT_DIR}/build $maus_coverage_extra \
+--exclude-dir build/test_maus_cpp/test_optics/
+
+python ${MAUS_ROOT_DIR}/src/map/MapCppGlobalTrackMatching/_test_MapCppGlobalTrackMatching.py
+
+echo 'Testing optics separately as it corrupts globals singleton'
+python ${MAUS_THIRD_PARTY}/third_party/install/bin/nosetests -v ${MAUS_ROOT_DIR}/build/test_maus_cpp/test_optics/ $maus_coverage_extra
+
 if [ $maus_lcov ]; then
     if [ $maus_lcov -ne "0" ]; then
         echo Building lcov output
@@ -50,8 +57,6 @@
     fi
 fi
 
-python ${MAUS_ROOT_DIR}/src/map/MapCppGlobalTrackMatching/_test_MapCppGlobalTrackMatching.py
-
 cd $here
 
 echo

=== modified file 'third_party/bash/21root.bash'
--- third_party/bash/21root.bash	2016-10-24 12:31:03 +0000
+++ third_party/bash/21root.bash	2017-06-08 17:10:31 +0000
@@ -59,11 +59,18 @@
       tar xvfz ${MAUS_ROOT_DIR}/third_party/source/${filename} -C ${MAUS_ROOT_DIR}/third_party/source > /dev/null
       mv ${MAUS_ROOT_DIR}/third_party/source/root ${MAUS_ROOT_DIR}/third_party/source/${directory}
 
+
+      echo
+      echo "INFO: Patching:"
+      echo
+      cd ${MAUS_ROOT_DIR}/third_party/source/${directory}
+      patch -p1 < ${MAUS_ROOT_DIR}/third_party/source/gcc6.patch
+
+
+      echo
+      echo "INFO: Configuring:"
+      echo
       cd ${MAUS_ROOT_DIR}/third_party/build/${directory}
-
-      echo
-      echo "INFO: Configuring:"
-      echo
       sleep 1
       # Tell cmake explicitly about third_party/install/include, third_party/install/lib and third_party/install/lib64  to look for headers and libs
       sed -i "s@include(MacroEnsureVersion)@include(MacroEnsureVersion)\ninclude_directories(${MAUS_ROOT_DIR}/third_party/install/include)\nlink_directories(${MAUS_ROOT_DIR}/third_party/install/lib ${MAUS_ROOT_DIR}/third_party/install/lib64)@" ${MAUS_ROOT_DIR}/third_party/source/${directory}/CMakeLists.txt

=== modified file 'third_party/bash/47cdb_cpp.bash'
--- third_party/bash/47cdb_cpp.bash	2017-04-25 18:17:08 +0000
+++ third_party/bash/47cdb_cpp.bash	2017-06-08 17:10:31 +0000
@@ -1,6 +1,6 @@
 #!/usr/bin/env bash
 
-version=1.0.3
+version=1.0.4
 directory="cdb-C++.v${version}"
 filename="cdb.client.api-cpp.v${version}.tgz"
 md5file=${filename}.md5
@@ -17,7 +17,7 @@
 
 # safer in any case to get the source from micemine
 # v1.0.3
-url=http://micewww.pp.rl.ac.uk/attachments/8738/${filename}
+url=https://micewww.pp.rl.ac.uk/attachments/download/8954/${filename}
 
 echo
 echo 'INFO: Installing third party library CDB-C++' $version

=== modified file 'third_party/build_all.bash'
--- third_party/build_all.bash	2017-04-23 22:27:43 +0000
+++ third_party/build_all.bash	2017-06-08 17:10:31 +0000
@@ -73,13 +73,14 @@
     # Make the numpy includes findable
     old_dir=`pwd`
     cd ${MAUS_THIRD_PARTY}/third_party/install/include
-    ln -s ../lib/python2.7/site-packages/numpy/core/include/numpy numpy
+    ln -sf ../lib/python2.7/site-packages/numpy/core/include/numpy numpy
     cd ${old_dir}
 
     ${MAUS_ROOT_DIR}/third_party/bash/51xboa.bash
     ${MAUS_ROOT_DIR}/third_party/bash/42libxml2.bash -j $MAUS_NUM_THREADS
     ${MAUS_ROOT_DIR}/third_party/bash/43libxslt.bash -j $MAUS_NUM_THREADS
     ${MAUS_ROOT_DIR}/third_party/bash/44cdb.bash
+    ${MAUS_ROOT_DIR}/third_party/bash/54flex.bash -j $MAUS_NUM_THREADS
     ${MAUS_ROOT_DIR}/third_party/bash/46gsoap.bash
     ${MAUS_ROOT_DIR}/third_party/bash/47cdb_cpp.bash
 
@@ -102,7 +103,6 @@
     # ${MAUS_ROOT_DIR}/third_party/bash/82heprep.bash
 
     # Doxygen for code documentation, requires a recent version of flex
-    ${MAUS_ROOT_DIR}/third_party/bash/54flex.bash -j $MAUS_NUM_THREADS
     ${MAUS_ROOT_DIR}/third_party/bash/55doxygen.bash -j $MAUS_NUM_THREADS
 
     # MAUS should now build okay - now for the test and execution environment

=== removed file 'third_party/source/cdb.client.api-cpp.v1.0.3.tgz.md5'
--- third_party/source/cdb.client.api-cpp.v1.0.3.tgz.md5	2017-04-25 15:09:46 +0000
+++ third_party/source/cdb.client.api-cpp.v1.0.3.tgz.md5	1970-01-01 00:00:00 +0000
@@ -1,1 +0,0 @@
-042024d29d62c23d1fa72f6002820574  cdb.client.api-cpp.v1.0.3.tgz

=== added file 'third_party/source/cdb.client.api-cpp.v1.0.4.tgz.md5'
--- third_party/source/cdb.client.api-cpp.v1.0.4.tgz.md5	1970-01-01 00:00:00 +0000
+++ third_party/source/cdb.client.api-cpp.v1.0.4.tgz.md5	2017-06-08 17:10:31 +0000
@@ -0,0 +1,1 @@
+6bd9f083ab1507eb41dd5c8e08ae0253  cdb.client.api-cpp.v1.0.4.tgz

=== added file 'third_party/source/gcc6.patch'
--- third_party/source/gcc6.patch	1970-01-01 00:00:00 +0000
+++ third_party/source/gcc6.patch	2017-06-08 17:10:31 +0000
@@ -0,0 +1,74 @@
+diff --git a/cint/cint/Module.mk b/cint/cint/Module.mk
+index 0293321..99a729b 100644
+--- a/cint/cint/Module.mk
++++ b/cint/cint/Module.mk
+@@ -165,6 +165,10 @@ ifeq ($(GCC_MAJOR),5)
+ CINTS2       := $(filter-out $(MODDIRSD)/libstrm.%,$(CINTS2))
+ CINTS2       += $(MODDIRSD)/gcc4strm.cxx
+ endif
++ifeq ($(GCC_MAJOR),6)
++CINTS2       := $(filter-out $(MODDIRSD)/libstrm.%,$(CINTS2))
++CINTS2       += $(MODDIRSD)/gcc4strm.cxx
++endif
+ ifneq ($(CLANG_MAJOR),)
+ CINTS2       := $(filter-out $(MODDIRSD)/libstrm.%,$(CINTS2))
+ CINTS2       += $(MODDIRSD)/gcc4strm.cxx
+@@ -207,6 +211,9 @@ IOSENUMC     := $(CINTDIRIOSEN)/iosenum.cxx
+ ifneq ($(CLANG_MAJOR),)
+ IOSENUMA     := $(CINTDIRIOSEN)/iosenum.$(ARCH)3
+ else
++ifeq ($(GCC_MAJOR),6)
++IOSENUMA     := $(CINTDIRIOSEN)/iosenum.$(ARCH)3
++else
+ ifeq ($(GCC_MAJOR),5)
+ IOSENUMA     := $(CINTDIRIOSEN)/iosenum.$(ARCH)3
+ else
+@@ -221,6 +228,7 @@ endif
+ endif
+ endif
+ endif
++endif
+ 
+ # used in the main Makefile
+ ALLHDRS     += $(CINTHT) $(CINTINCLUDES)
+@@ -355,6 +363,9 @@ endif
+ ifeq ($(GCC_MAJOR),5)
+ $(call stripsrc,$(CINTDIRSD)/gcc4strm.o): CINTCXXFLAGS += -Wno-strict-aliasing
+ endif
++ifeq ($(GCC_MAJOR),6)
++$(call stripsrc,$(CINTDIRSD)/gcc4strm.o): CINTCXXFLAGS += -Wno-strict-aliasing
++endif
+ 
+ 
+ $(MAKECINTO) $(CINTO): $(CINTCONF) $(ORDER_) $(CINTINCLUDES)
+@@ -389,8 +400,12 @@ endif
+ ##### configcint.h
+ ifeq ($(CPPPREP),)
+ # cannot use "CPPPREP?=", as someone might set "CPPPREP="
++ifeq ($(GCC_MAJOR),6)
++  CPPPREP = $(CXX) -std=c++98 -E -C
++else
+   CPPPREP = $(CXX) -E -C
+ endif
++endif
+ 
+ include $(CINTCONFMK)
+ ##### configcint.h - END
+diff --git a/cmake/modules/SetUpLinux.cmake b/cmake/modules/SetUpLinux.cmake
+index cc389b4..01b3b62 100644
+--- a/cmake/modules/SetUpLinux.cmake
++++ b/cmake/modules/SetUpLinux.cmake
+@@ -69,7 +69,12 @@ if(CMAKE_COMPILER_IS_GNUCXX)
+   set(CMAKE_C_FLAGS_PROFILE          "-g3 -fno-inline -ftest-coverage -fprofile-arcs")
+ 
+   #Settings for cint
+-  set(CPPPREP "${CXX} -E -C")
++  if (GCC_MAJOR EQUAL 6)
++    set(CPPPREP "${CXX} -std=c++98 -E -C")
++  else()
++    set(CPPPREP "${CXX} -E -C")
++  endif()
++
+   set(CXXOUT "-o ")
+   set(EXPLICITLINK "no") #TODO
+