diff --git a/UserTools/Factory/Factory.cpp b/UserTools/Factory/Factory.cpp index 7d1366c60..de800a3e3 100644 --- a/UserTools/Factory/Factory.cpp +++ b/UserTools/Factory/Factory.cpp @@ -176,5 +176,6 @@ if (tool=="AssignBunchTimingMC") ret=new AssignBunchTimingMC; if (tool=="FitRWMWaveform") ret=new FitRWMWaveform; if (tool=="PMTWaveformSim") ret=new PMTWaveformSim; if (tool=="LAPPDWaveformDisplay") ret=new LAPPDWaveformDisplay; +if (tool=="VertexLeastSquares") ret=new VertexLeastSquares; return ret; } diff --git a/UserTools/Unity.h b/UserTools/Unity.h index ea9fdd96d..7f72c5c30 100644 --- a/UserTools/Unity.h +++ b/UserTools/Unity.h @@ -184,3 +184,4 @@ #include "FitRWMWaveform.h" #include "PMTWaveformSim.h" #include "LAPPDWaveformDisplay.h" +#include "VertexLeastSquares.h" diff --git a/UserTools/VertexLeastSquares/MatrixUtil.cpp b/UserTools/VertexLeastSquares/MatrixUtil.cpp new file mode 100644 index 000000000..ff44c74e1 --- /dev/null +++ b/UserTools/VertexLeastSquares/MatrixUtil.cpp @@ -0,0 +1,209 @@ +#include "MatrixUtil.h" +#include + +namespace util +{ + //------------------------------------------------------------------------------ + void vector_show(const Vector& m, std::string str) + { + std::cout << std::fixed; + std::cout << std::setprecision(3); + for(int i = 0; i < m.size; i++) { + std::cout << m(i) << " "; + } + std::cout << "\n"; + } + + //------------------------------------------------------------------------------ + // c = a + b * s + void vmadd(const Vector& a, const Vector& b, double s, Vector& c) + { + if (c.size != a.size or c.size != b.size) { + std::cerr << "[vmadd]: vector sizes don't match\n"; + return; + } + + for (int i = 0; i < c.size; i++) + c(i) = a(i) + s * b(i); + } + + void mult_matvec(const Matrix& a, const Vector& b, Vector &c) + { + if (a.n != b.size) { + std::cerr << "[matvecmult]: Matrix-Vector multiplication not possible, sizes don't match !\n"; + return; + } + + if (a.m != c.size) { + std::cerr << "[matvecmult]: vector sizes don't match\n"; + return; + } + + for (int i = 0; i < a.m; i++) + for (int k = 0; k < a.n; k++) + c(i) += a(i,k) * b(k); + } + + //------------------------------------------------------------------------------ + // mat = I - 2*v*v^T + // !!! m is allocated here !!! + void compute_householder_factor(Matrix& mat, const Vector& v) + { + + int n = v.size; + mat.allocate(n,n); + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) + mat(i,j) = -2 * v(i) * v(j); + for (int i = 0; i < n; i++) + mat(i,i) += 1; + } + + //------------------------------------------------------------------------------ + // take c-th column of a matrix, put results in Vector v + void Matrix::extract_column(Vector& v, int c) { + if (m != v.size) { + std::cerr << "[Matrix::extract_column]: Matrix and Vector sizes don't match\n"; + return; + } + + for (int i = 0; i < m; i++) + v(i) = (*this)(i,c); + } + + //------------------------------------------------------------------------------ + void matrix_show(const Matrix& m, std::string str) + { + std::cout << std::fixed; + std::cout << std::setprecision(3); + for(int i = 0; i < m.m; i++) { + for (int j = 0; j < m.n; j++) { + std::cout << m(i,j) << " "; + } + std::cout << "\n"; + } + std::cout << "\n"; + } + + //------------------------------------------------------------------------------ + // L2-norm ||A-B||^2 + double matrix_compare(const Matrix& A, const Matrix& B) + { + // matrices must have same size + if (A.m != B.m or A.n != B.n) + return std::numeric_limits::max(); + + double res=0; + for(int i = 0; i < A.m; i++) { + for (int j = 0; j < A.n; j++) { + res += (A(i,j)-B(i,j)) * (A(i,j)-B(i,j)); + } + } + + res /= A.m*A.n; + return res; + } + + //------------------------------------------------------------------------------ + // Householder transformation to obtain QR matrices + void householder(const Matrix &mat, Matrix& R, Matrix& Q) + { + + int m = mat.m; + int n = mat.n; + + // array of factor Q1, Q2, ... Qm + std::vector qv(m); + + // temp array + Matrix z(mat); + Matrix z1; + + for (int k = 0; k < n && k < m - 1; k++) { + + Vector e(m), x(m); + double a; + + // compute minor + z1.compute_minor(z, k); + + // extract k-th column into x + z1.extract_column(x, k); + + a = x.norm(); + if (mat(k,k) > 0) a = -a; + + for (int i = 0; i < e.size; i++) + e(i) = (i == k) ? 1 : 0; + + // e = x + a*e + vmadd(x, e, a, e); + + // e = e / ||e|| + e.rescale_unit(); + + // qv[k] = I - 2 *e*e^T + compute_householder_factor(qv[k], e); + + // z = qv[k] * z1 + z.mult(qv[k], z1); + + } + + Q = qv[0]; + + // after this loop, we will obtain Q (up to a transpose operation) + for (int i = 1; i < n && i < m - 1; i++) { + + z1.mult(qv[i], Q); + Q = z1; + + } + + R.mult(Q, mat); + Q.transpose(); + } + + //------------------------------------------------------------------------------ + void solve_upper_triangular(const Matrix &r, const Vector &b, Vector &solution) + { + const int32_t ncolumns = solution.size; + + for ( int32_t k = ncolumns - 1; k >= 0; --k ) { + double total = 0.0; + for ( int32_t j = k + 1; j < ncolumns; ++j ) { + total += r(k, j) * solution(j); + } + solution(k) =( b(k) - total ) / r(k, k); + } + } + + //------------------------------------------------------------------------------ + void least_squares(const Matrix &A, const Vector &b, Vector& solution, bool verbose) + { + Matrix Q(A.n, A.n); + Matrix R(A.m, A.n); + + householder(A, R, Q); + + Q.transpose(); + + Vector c(Q.m); + mult_matvec(Q, b, c); + + if (verbose) { + std::cout << "A" << std::endl; + matrix_show(A, ""); + std::cout << "b" << std::endl; + vector_show(b, ""); + std::cout << "Q" << std::endl; + matrix_show(Q, ""); + std::cout << "R" << std::endl; + matrix_show(R, ""); + std::cout << "c" << std::endl; + vector_show(c, ""); + } + + solve_upper_triangular(R, c, solution); + } +} diff --git a/UserTools/VertexLeastSquares/MatrixUtil.h b/UserTools/VertexLeastSquares/MatrixUtil.h new file mode 100644 index 000000000..bd91de7c0 --- /dev/null +++ b/UserTools/VertexLeastSquares/MatrixUtil.h @@ -0,0 +1,254 @@ +#include +#include +#include // for memset +#include +#include +#include + +#include + +namespace util +{ + + // column vector + class Vector { + + public: + // default constructor (don't allocate) + Vector() : size(0), data(nullptr) {} + + // constructor with memory allocation, initialized to zero + Vector(int size_) : Vector() { + size = size_; + allocate(size_); + } + + // destructor + ~Vector() { + deallocate(); + } + + // access data operators + double& operator() (int i) { + return data[i]; } + double operator() (int i) const { + return data[i]; } + + // operator assignment + Vector& operator=(const Vector& source) { + + // self-assignment check + if (this != &source) { + if ( size != (source.size) ) { // storage cannot be reused + allocate(source.size); // re-allocate storage + } + // storage can be used, copy data + std::copy(source.data, source.data + source.size, data); + } + return *this; + } + + // memory allocation + void allocate(int size_) { + + deallocate(); + + // new sizes + size = size_; + + data = new double[size_]; + memset(data,0,size_*sizeof(double)); + + } // allocate + + // memory free + void deallocate() { + + if (data) + delete[] data; + + data = nullptr; + + } + + // ||x|| + double norm() { + double sum = 0; + for (int i = 0; i < size; i++) sum += (*this)(i) * (*this)(i); + return sqrt(sum); + } + + // divide data by factor + void rescale(double factor) { + for (int i = 0; i < size; i++) (*this)(i) /= factor; + } + + void rescale_unit() { + double factor = norm(); + rescale(factor); + } + + int size; + + private: + double* data; + + }; // class Vector + + //------------------------------------------------------------------------------ + class Matrix { + + public: + // default constructor (don't allocate) + Matrix() : m(0), n(0), data(nullptr) {} + + // constructor with memory allocation, initialized to zero + Matrix(int m_, int n_) : Matrix() { + m = m_; + n = n_; + allocate(m_,n_); + } + + // copy constructor + Matrix(const Matrix& mat) : Matrix(mat.m,mat.n) { + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + (*this)(i,j) = mat(i,j); + } + + // constructor from array + template + Matrix(double (&a)[rows][cols]) : Matrix(rows,cols) { + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + (*this)(i,j) = a[i][j]; + } + + // constructor from Vector, creates a vertical matrix + Matrix(const Vector& vec) : Matrix(vec.size, 1) { + + for (int i = 0; i < m; i++) + (*this)(i,0) = vec(i); + } + + + // destructor + ~Matrix() { + deallocate(); + } + + + // access data operators + double& operator() (int i, int j) { + return data[i+m*j]; } + double operator() (int i, int j) const { + return data[i+m*j]; } + + // operator assignment + Matrix& operator=(const Matrix& source) { + + // self-assignment check + if (this != &source) { + if ( (m*n) != (source.m * source.n) ) { // storage cannot be reused + allocate(source.m,source.n); // re-allocate storage + } + // storage can be used, copy data + std::copy(source.data, source.data + source.m*source.n, data); + } + return *this; + } + + // compute minor + void compute_minor(const Matrix& mat, int d) { + + allocate(mat.m, mat.n); + + for (int i = 0; i < d; i++) + (*this)(i,i) = 1.0; + for (int i = d; i < mat.m; i++) + for (int j = d; j < mat.n; j++) + (*this)(i,j) = mat(i,j); + + } + + // Matrix multiplication + // c = a * b + // c will be re-allocated here + void mult(const Matrix& a, const Matrix& b) { + + if (a.n != b.m) { + std::cerr << "Matrix multiplication not possible, sizes don't match !\n"; + return; + } + + // reallocate ourself if necessary i.e. current Matrix has not valid sizes + if (a.m != m or b.n != n) + allocate(a.m, b.n); + + memset(data,0,m*n*sizeof(double)); + + for (int i = 0; i < a.m; i++) + for (int j = 0; j < b.n; j++) + for (int k = 0; k < a.n; k++) + (*this)(i,j) += a(i,k) * b(k,j); + + } + + void transpose() { + for (int i = 0; i < m; i++) { + for (int j = 0; j < i; j++) { + double t = (*this)(i,j); + (*this)(i,j) = (*this)(j,i); + (*this)(j,i) = t; + } + } + } + + // take c-th column of m, put in v + void extract_column(Vector& v, int c); + + // memory allocation + void allocate(int m_, int n_) { + + // if already allocated, memory is freed + deallocate(); + + // new sizes + m = m_; + n = n_; + + data = new double[m_*n_]; + memset(data,0,m_*n_*sizeof(double)); + + } // allocate + + // memory free + void deallocate() { + + if (data) + delete[] data; + + data = nullptr; + + } + + int m, n; + + private: + double* data; + + }; // struct Matrix + + // Other useful functions + void vector_show(const Vector& m, std::string str); + void vmadd(const Vector& a, const Vector& b, double s, Vector& c); + void mult_matvec(const Matrix& a, const Vector& b, Vector &c); + void compute_householder_factor(Matrix& mat, const Vector& v); + void matrix_show(const Matrix& m, std::string str); + double matrix_compare(const Matrix& A, const Matrix& B); + void householder(const Matrix &mat, Matrix& R, Matrix& Q); + void solve_upper_triangular(const Matrix &r, const Vector &b, Vector &solution); + void least_squares(const Matrix &A, const Vector &b, Vector &solution, bool verbose = false); +} diff --git a/UserTools/VertexLeastSquares/README.md b/UserTools/VertexLeastSquares/README.md new file mode 100644 index 000000000..a351ac40f --- /dev/null +++ b/UserTools/VertexLeastSquares/README.md @@ -0,0 +1,102 @@ +# VertexLeastSquares + +VertexLeastSquares is a vertex reconstruction tool that relies only on PMT hit timing. It recommended for use in low energy events that do not have an MRD track to assist in reconstruction. It utilizes a matrix least squares method for vertexing interactions. It also does hit filtering / cleaning prior to reconstruction to only fit the set of PMT hits compatible with the assumption all light in the event originated from a single point. Prior to fitting, it populates the tank with seeds, then uses the Gauss-Newton method to update an initial guess until a global minimum is found. The tool will return the best-fit vertex position to the available PMT hits information. + +## Data + +**VertexLeastSquaresMap** `std::map` +* Stored in `ANNIEEvent`, the reconstructed vertex (x,y,z) [m] for a cluster. The vertices share the same coordinate system as WCSim. +* To access the reconstructed vertices: +``` +# grab vertex map +m_data->Stores.at("ANNIEEvent")->Get("VertexLeastSquaresMap", fVertexMap); + +# index the map with the cluster time for a given cluster +Position vertex = fVertexMap->at(cluster_time); # Position is a ToolAnalysis class + +# grab vertex +recoVtxX = vertex.X() // [m] +recoVtxY = vertex.Y() // [m] +recoVtxZ = vertex.Z() // [m] +``` +* To shift the vertex positions back to a coordinate system with the ANNIE tank center at (0,0,0): +``` +recoVtxX = (vertex.X() + 0) // [m] +recoVtxY = (vertex.Y() + 0.1446) // [m] +recoVtxZ = (vertex.Z() - 1.681) // [m] +``` + +Note that the reconstructed vertex is given as the best fit to the hit times within a cluster. When comparing this information to GENIE, any reconstructed, delayed clusters will necessarily be "off" from the true interaction vertex. If neutrons, the reconstructed vertex could be compared to the neutron capture point instead. + +If a suitable vertex is not found, the map is populated with `-99`. + +**VertexLeastSquaresStdevMap** `std::map` +* Stored in `ANNIEEvent`, the stdev from the timing residuals of the fitted vertex [unitless]. + +* To access the stdev of the fitted vertex: +``` +# grab stdev +m_data->Stores.at("ANNIEEvent")->Get("VertexLeastSquaresStdevMap", fVertexStdevMap); + +# index the map with the cluster time for a given cluster +fVertexStdevMap->at(cluster_time); +``` + +The stdev will also get a value of `-99` if a suitable vertex is not found. + +## Configuration + +``` +verbosity 0 +UseMCHits 0 # MC hits (1) or data (0). For PMTWaveformSim, use data (0) + +BreakDist 0.05 # how close two guesses have to be in order to end the vertex search [m] +Regularizer 1 # Tikhonov (L2) regularization term --> controls how far vertices are allowed to wander away from the center +FittingSteps 100 # how many iterations the algorithm will run through + +YSpacing 0.5 # vertical spacing of the sunflower seed layers [m] +NPlanarPoints 20 # how many seeds are arranged in a sunflower pattern per layer + +FilterHitsCompatibleWindow 0 # (0) = keep first 4 hits as "seed" hits, then keep later hits only if they are casually disconnected + # (1) = search for the most compatible (causally disconnected) set of hits with a sliding time window + +OnlyUseReliablePMTs 0 # mask out PMTs with poor timing uncertainties (determined from the laser analysis) +MinimumReliableTimingStd 1.5 # maximum timing uncertainty for PMT masking + +ExternalSeeding 0 # Add extend seeding outside tank (to help reconstruct external events), and allow guesses to wander outside tank volume +YBuffer 1.0 # additional (Y) buffer height added to the seeding volume (top and bottom, with respect to center) [m] +RadialBuffer 1.0 # additional radial distance added to the seeding volume (with respect to center) [m] +ExternalVertexMaxRadius 4.0 # maximum allowable radial distance a guess is allowed to be from the tank center (if external seeding is enabled) [m] +ExternalVertexMaxHeight 4.0 # maximum allowable Y distance a guess is allowed to be from the tank center (if external seeding is enabled) [m] +``` + +## Working toolchain + +For the parametric MC hits: +``` +LoadGeometry +LoadWCSim +ClusterFinder +VertexLeastSquares +``` + +For the waveform MC hits: +``` +LoadGeometry +LoadWCSim +PMTWaveformSim +PhaseIIADCHitFinder +ClusterFinder +VertexLeastSquares +``` + +For the data: +``` +LoadANNIEEvent +LoadGeometry +ClusterFinder +VertexLeastSquares +``` + +### Additional Information +More information on the tool can be found here: https://annie-docdb.fnal.gov/cgi-bin/sso/ShowDocument?docid=5621 \ No newline at end of file diff --git a/UserTools/VertexLeastSquares/VertexLeastSquares.cpp b/UserTools/VertexLeastSquares/VertexLeastSquares.cpp new file mode 100644 index 000000000..06c45901c --- /dev/null +++ b/UserTools/VertexLeastSquares/VertexLeastSquares.cpp @@ -0,0 +1,888 @@ +#include "VertexLeastSquares.h" + +VertexLeastSquares::VertexLeastSquares():Tool(){} + +//------------------------------------------------------------------------------ +bool VertexLeastSquares::Initialise(std::string configfile, DataModel &data) +{ + + /////////////////// Useful header /////////////////////// + if(configfile!="") m_variables.Initialise(configfile); // loading config file + //m_variables.Print(); + + m_data= &data; //assigning transient data pointer + ///////////////////////////////////////////////////////////////// + + // Load my config parameters + bool gotVerbosity = m_variables.Get("verbosity",verbosity); + if (!gotVerbosity) { + verbosity = 0; + Log("VertexLeastSquares: \"verbosity\" not set in the config, defaulting to 0", v_error, verbosity); + } + + // *********************************************************************** + + bool gotUseMCHits = m_variables.Get("UseMCHits", fUseMCHits); + if (!gotUseMCHits) { + Log("VertexLeastSquares: \"UseMCHits\" not set in the config! Aborting!", v_error, verbosity); + return false; + } + + bool gotBreakDist = m_variables.Get("BreakDist", fBreakDist); + if (!gotBreakDist) { + fBreakDist = 0.05; + Log("VertexLeastSquares: \"BreakDist\" not set in the config! Using default value of 5 cm.", v_error, verbosity); + } + + bool gotYSpacing = m_variables.Get("YSpacing", fYSpacing); + if (!gotYSpacing) { + fYSpacing = 0.5; + Log("VertexLeastSquares: \"YSpacing\" not set in the config! Using default value of 50 cm.", v_error, verbosity); + } + + bool gotNPlanarPoints = m_variables.Get("NPlanarPoints", fNPlanarPoints); + if (!gotNPlanarPoints) { + fNPlanarPoints = 20; + Log("VertexLeastSquares: \"NPlanarPoints\" not set in the config! Using default value of 20.", v_error, verbosity); + } + + bool gotRegularizer = m_variables.Get("Regularizer", fRegularizer); + if (!gotRegularizer) { + fRegularizer = 1.; + Log("VertexLeastSquares: \"Regularizer\" not set in the config! Using default value of 1.", v_error, verbosity); + } + + bool gotFittingSteps = m_variables.Get("FittingSteps", fFittingSteps); + if (!gotFittingSteps) { + fFittingSteps = 100.; + Log("VertexLeastSquares: \"FittingSteps\" not set in the config! Using default value of 100.", v_error, verbosity); + } + + bool gotMostCompatible = m_variables.Get("FilterHitsCompatibleWindow", fCompatibleHits); + if (!gotMostCompatible) { + fCompatibleHits = false; + Log("VertexLeastSquares: \"FilterHitsCompatibleWindow\" not set in the config! Defaulting to false.", v_error, verbosity); + } + + bool gotGoodTimes = m_variables.Get("OnlyUseReliablePMTs", fGoodTimes); + if (!gotGoodTimes) { + fGoodTimes = false; + Log("VertexLeastSquares: \"OnlyUseReliablePMTs\" not set in the config! Defaulting to false.", v_error, verbosity); + } + + bool gotMinGoodTime = m_variables.Get("MinimumReliableTimingStd", fMinGoodTime); + if (!gotMinGoodTime) { + fMinGoodTime = 1.5; // ns + Log("VertexLeastSquares: \"MinimumReliableTimingStd\" not set in the config! Defaulting to 1.5ns.", v_error, verbosity); + } + + // *********************************************************************** + + bool gotExternalSeeding = m_variables.Get("ExternalSeeding",fExternalSeeding); + if (!gotExternalSeeding) { + fExternalSeeding = false; + Log("VertexLeastSquares: \"ExternalSeeding\" not set in the config, defaulting to false (only seed the tank).", v_error, verbosity); + } + + bool gotYbuffer = m_variables.Get("YBuffer", fYBuffer); + if (!gotYbuffer) { + fYSpacing = 1.0; + Log("VertexLeastSquares: \"YBuffer\" not set in the config! Using default value of 100 cm.", v_error, verbosity); + } + + bool gotRadialbuffer = m_variables.Get("RadialBuffer", fRadialBuffer); + if (!gotRadialbuffer) { + fRadialBuffer = 1.0; + Log("VertexLeastSquares: \"RadialBuffer\" not set in the config! Using default value of 100 cm.", v_error, verbosity); + } + + bool gotVertexMaxRadius = m_variables.Get("ExternalVertexMaxRadius", fExternalVertexMaxRadius); + if (!gotVertexMaxRadius) { + fExternalVertexMaxRadius = 4.0; + Log("VertexLeastSquares: \"ExternalVertexMaxRadius\" not set in the config! Using default value of 4 m.", v_error, verbosity); + } + + bool gotVertexMaxHeight = m_variables.Get("ExternalVertexMaxHeight", fExternalVertexMaxHeight); + if (!gotVertexMaxHeight) { + fExternalVertexMaxHeight = 4.0; + Log("VertexLeastSquares: \"ExternalVertexMaxHeight\" not set in the config! Using default value of 4 m.", v_error, verbosity); + } + + // *********************************************************************** + + fNBoundary = int(sqrt(fRegularizer)); + + // Load the geometry service + bool gotGeometry = m_data->Stores.at("ANNIEEvent")->Header->Get("AnnieGeometry",fGeom); + if(!gotGeometry){ + Log("VertexLeastSquares: Error retrieving Geometry from ANNIEEvent! Aborting!", v_error, verbosity); + return false; + } + + // Load the PMT timing uncertainties + if (fGoodTimes) { + bool got_timing_std = m_data->CStore.Get("ChannelNumToTankPMTTimingSigmaMap",ChannelKeyToTimingSigmaMap); + if (!got_timing_std) { + logmessage = "VertexLeastSquares: \"OnlyUseReliablePMTs\" enabled --> error retrieving PMT timing uncertainty... double check LoadGeometry maybe?"; + Log(logmessage, v_error, verbosity); + return false; + } + } + + // Set up the pointers we're going to save. No need to + // delete it at Finalize, the store will handle it (I'm trusting you Andrew!) + fVertexMap = new std::map; + fVertexStdevMap = new std::map; + + return true; +} + +//------------------------------------------------------------------------------ +bool VertexLeastSquares::Execute() +{ + + fVertexMap->clear(); + fVertexStdevMap->clear(); + if (fUseMCHits) { + bool gotClusters = m_data->CStore.Get("ClusterMapMC", fClusterMapMC); + if (!gotClusters) { + logmessage = "VertexLeastSquares: no ClusterMapMC in the ANNIEEvent!!! Are you sure you ran ClusterFinder?"; + Log(logmessage, v_error, verbosity); + return false; + } + RunLoopMC(); + } else { + bool gotClusters = m_data->CStore.Get("ClusterMap", fClusterMap); + if (!gotClusters) { + logmessage = "VertexLeastSquares: no ClusterMap in the ANNIEEvent!!! Are you sure you ran ClusterFinder?"; + Log(logmessage, v_error, verbosity); + return false; + } + + RunLoop(); + } + + m_data->Stores.at("ANNIEEvent")->Set("VertexLeastSquaresMap", fVertexMap); + m_data->Stores.at("ANNIEEvent")->Set("VertexLeastSquaresStdevMap", fVertexStdevMap); + + return true; +} + +//------------------------------------------------------------------------------ +bool VertexLeastSquares::Finalise() +{ + + return true; +} + +//----------------------------------------------------------------------------- +std::vector VertexLeastSquares::GenerateVetices() // this is Andrew's fault there's a typo (and no I'm not going to fix it :) ) +{ + + if (fExternalSeeding) { + + // generate vertices in and outside the tank (+ some buffer) + + std::vector vertices; + + std::vector ys; + double yMin = fGeom->GetTankCentre().Y() - fGeom->GetTankHalfheight() - fYBuffer; // tank height +/- Y buffer assigned in config + double yMax = fGeom->GetTankCentre().Y() + fGeom->GetTankHalfheight() + fYBuffer; + for (double yValue = yMin; yValue <= yMax; yValue += fYSpacing) + ys.push_back(yValue); + + // generate the x and z values using the sunflower pattern (extending to the radius + buffer) + std::vector xs; + std::vector zs; + double max_radius = fGeom->GetTankRadius() + fRadialBuffer; + for (int n = 1; n < fNPlanarPoints+1; ++n) { + double rad = ( (n > fNPlanarPoints + fNBoundary) ? 1.0 : + sqrt((n+0.5)/(fNPlanarPoints - (fNBoundary+1.)/2.)) ); + + // scale to the extended radius + rad *= max_radius; + rad = rad > max_radius ? max_radius : rad; + + double angle = 2. * pi * n / phisq; + + xs.push_back(rad*cos(angle)); + zs.push_back(rad*sin(angle)); + } + + // add vertices from the sunflower layers + for (uint yIdx = 0; yIdx < ys.size(); ++yIdx) { + + // Apply a random rotation to each y level to cover more space + double rotAng = (double(rand())/RAND_MAX)*2*pi; + for (uint xzIdx = 0; xzIdx < fNPlanarPoints; ++xzIdx) { + double x = xs[xzIdx]; + double z = zs[xzIdx]; + + vertices.push_back(Position(x*cos(rotAng) - z*sin(rotAng) + fGeom->GetTankCentre().X(), + ys[yIdx], + z*cos(rotAng) + x*sin(rotAng) + fGeom->GetTankCentre().Z())); + } + } + + return vertices; + + } else { + + // generate vertices within the tank + std::vector ys; + double yMin = fGeom->GetTankCentre().Y() - fGeom->GetTankHalfheight(); + double yMax = fGeom->GetTankCentre().Y() + fGeom->GetTankHalfheight(); + for (double yValue = yMin; yValue <= yMax; yValue += fYSpacing) + ys.push_back(yValue); + + // generate the x and z values using the sunflower pattern + std::vector xs; + std::vector zs; + for (int n = 1; n < fNPlanarPoints+1; ++n) { + double rad = ( (n > fNPlanarPoints + fNBoundary) ? 1.0 : + sqrt((n+0.5)/(fNPlanarPoints - (fNBoundary+1.)/2.)) ); + + // Scale the radius to the tank and make sure it's inside + rad *= fGeom->GetTankRadius(); + rad = rad > fGeom->GetTankRadius() ? fGeom->GetTankRadius() : rad; + + double angle = 2. * pi * n / phisq; + + xs.push_back(rad*cos(angle)); + zs.push_back(rad*sin(angle)); + } + + std::vector vertices; + for (uint yIdx = 0; yIdx < ys.size(); ++yIdx) { + + // Apply a random rotation to each y level to cover more space + double rotAng = (double(rand())/RAND_MAX)*2*pi; + for (uint xzIdx = 0; xzIdx < fNPlanarPoints; ++xzIdx) { + double x = xs[xzIdx]; + double z = zs[xzIdx]; + + vertices.push_back(Position(x*cos(rotAng) - z*sin(rotAng) + fGeom->GetTankCentre().X(), + ys[yIdx], + z*cos(rotAng) + x*sin(rotAng) + fGeom->GetTankCentre().Z())); + } + } + + return vertices; + + } +} + +//------------------------------------------------------------------------------ +void VertexLeastSquares::EvalAtGuessVertex(util::Matrix &A, util::Vector &b, + const Position &guess, + const std::vector &hits) +{ + // The function f(x) = 0 + // 0 = 1/c * sqrt( (X-xi)^2 + (Y-yi)^2 + (Z-zi)^2 ) + T - ti + // A = Jacobian(guess), b = -f(guess) + double mean_T = 0; + for (uint hitIdx = 0; hitIdx < hits.size(); ++hitIdx) { + Detector *det = fGeom->ChannelToDetector(hits[hitIdx].GetTubeId()); + Position detpos = det->GetDetectorPosition(); + double dist = (detpos - guess).Mag(); + double ti = hits[hitIdx].GetTime(); + + A(hitIdx, 0) = (guess.X() - detpos.X())/(dist * fSoL); + A(hitIdx, 1) = (guess.Y() - detpos.Y())/(dist * fSoL); + A(hitIdx, 2) = (guess.Z() - detpos.Z())/(dist * fSoL); + A(hitIdx, 3) = 1.; + + b(hitIdx) = -dist/fSoL + ti; + + mean_T += ti - dist/fSoL; + }// end loop over hits + + mean_T = mean_T / hits.size(); + + // Subtract off the mean emission time (T) + for (uint hitIdx = 0; hitIdx < hits.size(); ++hitIdx) + b(hitIdx) -= mean_T; + + // Tack on the regularization term + A(hits.size(), 0) = fRegularizer; + A(hits.size() + 1, 1) = fRegularizer; + A(hits.size() + 2, 2) = fRegularizer; + A(hits.size() + 3, 3) = fRegularizer; + b(hits.size() ) = 0; + b(hits.size() + 1) = 0; + b(hits.size() + 2) = 0; + b(hits.size() + 3) = 0; + + +} + + +//------------------------------------------------------------------------------ +void VertexLeastSquares::EvalAtGuessVertexMC(util::Matrix &A, util::Vector &b, + const Position &guess, + const std::vector &hits) +{ + // The function f(x) = 0 + // 0 = 1/c * sqrt( (X-xi)^2 + (Y-yi)^2 + (Z-zi)^2 ) + T - ti + // A = Jacobian(guess), b = -f(guess) + double mean_T = 0; + for (uint hitIdx = 0; hitIdx < hits.size(); ++hitIdx) { + Detector *det = fGeom->ChannelToDetector(hits[hitIdx].GetTubeId()); + Position detpos = det->GetDetectorPosition(); + double dist = (detpos - guess).Mag(); + double ti = hits[hitIdx].GetTime(); + + A(hitIdx, 0) = (guess.X() - detpos.X())/(dist * fSoL); + A(hitIdx, 1) = (guess.Y() - detpos.Y())/(dist * fSoL); + A(hitIdx, 2) = (guess.Z() - detpos.Z())/(dist * fSoL); + A(hitIdx, 3) = 1.; + + b(hitIdx) = -dist/fSoL + ti; + + mean_T += ti - dist/fSoL; + }// end loop over hits + + mean_T = mean_T / hits.size(); + + // Subtract off the mean emission time (T) + for (uint hitIdx = 0; hitIdx < hits.size(); ++hitIdx) + b(hitIdx) -= mean_T; + + // Tack on the regularization term + A(hits.size(), 0) = fRegularizer; + A(hits.size() + 1, 1) = fRegularizer; + A(hits.size() + 2, 2) = fRegularizer; + A(hits.size() + 3, 3) = fRegularizer; + b(hits.size() ) = 0; + b(hits.size() + 1) = 0; + b(hits.size() + 2) = 0; + b(hits.size() + 3) = 0; + + +} + +//------------------------------------------------------------------------------ +void VertexLeastSquares::RunLoop() +{ + + for (auto clusterpair : *fClusterMap) { + // Filter hits + auto filt_hits = FilterHits(clusterpair.second); + + // We need at least 4 hits + if (filt_hits.size() < 4) { + std::cout << "Not enough hits. We only have: " << filt_hits.size() << std::endl; + fVertexMap->emplace(clusterpair.first, Position(-99, -99, -99)); + fVertexStdevMap->emplace(clusterpair.first, -99); + continue; + } + + // Loop over vertex seeds + std::vector seeds = GenerateVetices(); + + std::vector bestVtxs; + double bestStDev = 999999; + int bestIdx = 0; + int seednum = 0; + for (auto guess : seeds) { + Position lastGuess; + int loop_num = 0; + bool inTank = true; + auto hits = filt_hits; + while (true) { + lastGuess = guess; + inTank = true; + + // Set up Matrix and solution vector + util::Matrix A(filt_hits.size() + 4, 4); + util::Vector b(filt_hits.size() + 4); + EvalAtGuessVertex(A, b, guess, hits); + + // Solve it and update the guess vertex + util::Vector solution(4); + util::least_squares(A, b, solution); + + // If the solution has nans it's because all rows of A evaluate to the same thing + // revert to the last vertex and break + if (std::isnan(solution(0)) || std::isnan(solution(1)) || std::isnan(solution(2))) { + inTank = false; + break; + + } + + // Otherwise, update the guess and carry on + guess = lastGuess + Position(solution(0), solution(1), solution(2)); + + // Exit if the new guess is close enough to the last one + if ((guess - lastGuess).Mag() < fBreakDist) break; + + if (!fExternalSeeding) { + // Exit if we're outside of the tank + if(!fGeom->GetTankContained(guess)) { + inTank = false; + break; + } + } + + // Don't let the loop last forever + ++loop_num; + if (loop_num > fFittingSteps) break; + }// end while loop + + // even tho we are seeding guesses outside the tank, still check if the guess is "reasonable" + // the tank-only setting just calls GetTankContained(guess) - we instead must check if the guess is within our defined max radius and half height + // "reasonable" distance away for a guess: radius, height +/- ExternalVertexMax(Radius)Height defined in config + if (fExternalSeeding) { + + Position tankCenter = fGeom->GetTankCentre(); + double dx = guess.X() - tankCenter.X(); + double dz = guess.Z() - tankCenter.Z(); + double dy = guess.Y() - tankCenter.Y(); + double radial_dist = std::sqrt(dx*dx + dz*dz); + double max_radial = fGeom->GetTankRadius() + fExternalVertexMaxRadius; + double max_vertical = fGeom->GetTankHalfheight() + fExternalVertexMaxHeight; + + if (radial_dist > max_radial || std::abs(dy) > max_vertical) { + continue; // reject this seed + } + + } + + // Skip if no tank contained vertex was found for this seed + if (!inTank) continue; + + bestVtxs.push_back(guess); + + // Check the std deviation of b for the best guess from this seed + util::Matrix A(filt_hits.size() + 4, 4); + util::Vector b(filt_hits.size() + 4); + EvalAtGuessVertex(A, b, bestVtxs.back(), filt_hits); + double sum = 0; + double sum_sq = 0; + for (uint idx = 0; idx < b.size - 4; ++idx) { + sum += b(idx); + sum_sq += b(idx) * b(idx) ; + } + double mean = sum / b.size; + double stdev = sum_sq / b.size - mean * mean; + if (stdev < bestStDev) { + bestStDev = stdev; + bestIdx = bestVtxs.size()-1; + } + ++seednum; + }// end loop over seed vertices + + // Save the vertex and stdev if we have one, otherwise default to -99s + if (bestVtxs.size()) { + fVertexMap->emplace(clusterpair.first, bestVtxs[bestIdx]); + fVertexStdevMap->emplace(clusterpair.first, bestStDev); + + } else { + fVertexMap->emplace(clusterpair.first, Position(-99, -99, -99)); + fVertexStdevMap->emplace(clusterpair.first, -99); + } + + }// end loop over the clusters +} + +//------------------------------------------------------------------------------ +void VertexLeastSquares::RunLoopMC() +{ + for (auto clusterpair : *fClusterMapMC) { + // Filter hits + auto filt_hits = FilterHitsMC(clusterpair.second); + + //int PMTID = mcHitsIt.first; + + // We need at least 4 hits + if (filt_hits.size() < 4) { + std::cout << "Not enough hits. We only have: " << filt_hits.size() << std::endl; + fVertexMap->emplace(clusterpair.first, Position(-99, -99, -99)); + fVertexStdevMap->emplace(clusterpair.first, -99); + continue; + } + + // Loop over vertex seeds + std::vector seeds = GenerateVetices(); + + std::vector bestVtxs; + double bestStDev = 999999; + int bestIdx = 0; + int seednum = 0; + for (auto guess : seeds) { + Position lastGuess; + int loop_num = 0; + bool inTank = true; + auto hits = filt_hits; + while (true) { + lastGuess = guess; + inTank = true; + + // Set up Matrix and solution vector + util::Matrix A(hits.size() + 4, 4); + util::Vector b(hits.size() + 4); + EvalAtGuessVertexMC(A, b, guess, hits); + + // Solve it and update the guess vertex + util::Vector solution(4); + util::least_squares(A, b, solution); + + // If the solution has nans it's because all rows of A evaluate to the same thing + // revert to the last vertex and break + if (std::isnan(solution(0)) || std::isnan(solution(1)) || std::isnan(solution(2))) { + inTank = false; + break; + + } + + // Otherwise, update the guess and carry on + guess = lastGuess + Position(solution(0), solution(1), solution(2)); + + // Exit if the new guess is close enough to the last one + if ((guess - lastGuess).Mag() < fBreakDist) break; + + if (!fExternalSeeding) { + // Exit if we're outside of the tank + if(!fGeom->GetTankContained(guess)) { + inTank = false; + break; + } + } + + // Don't let the loop last forever + ++loop_num; + if (loop_num > fFittingSteps) break; + }// end while loop + + // even tho we are seeding guesses outside the tank, still check if the guess is "reasonable" + // the tank-only setting just calls GetTankContained(guess) - we instead must check if the guess is within our defined max radius and half height + // "reasonable" distance away for a guess: radius, height +/- ExternalVertexMax(Radius)Height defined in config + if (fExternalSeeding) { + + Position tankCenter = fGeom->GetTankCentre(); + double dx = guess.X() - tankCenter.X(); + double dz = guess.Z() - tankCenter.Z(); + double dy = guess.Y() - tankCenter.Y(); + double radial_dist = std::sqrt(dx*dx + dz*dz); + double max_radial = fGeom->GetTankRadius() + fExternalVertexMaxRadius; + double max_vertical = fGeom->GetTankHalfheight() + fExternalVertexMaxHeight; + + if (radial_dist > max_radial || std::abs(dy) > max_vertical) { + continue; // reject this seed + } + + } + + // Skip if no tank contained vertex was found for this seed + if (!inTank) continue; + + // Check the std deviation of b for the best guess from this seed + bestVtxs.push_back(guess); + util::Matrix A(filt_hits.size() + 4, 4); + util::Vector b(filt_hits.size() + 4); + EvalAtGuessVertexMC(A, b, bestVtxs.back(), filt_hits); + double sum = 0; + double sum_sq = 0; + for (uint idx = 0; idx < b.size - 4; ++idx) { + sum += b(idx); + sum_sq += b(idx) * b(idx) ; + } + double mean = sum / b.size; + double stdev = sum_sq / b.size - mean * mean; + if (stdev < bestStDev) { + bestStDev = stdev; + bestIdx = bestVtxs.size()-1; + } + ++seednum; + }// end loop over seed vertices + + // Save the vertex if we have one, otherwise default to -99s + if (bestVtxs.size()) { + fVertexMap->emplace(clusterpair.first, bestVtxs[bestIdx]); + fVertexStdevMap->emplace(clusterpair.first, -99); + + } + else { + std::cout << "No good vertex found!!" << std::endl; + fVertexMap->emplace(clusterpair.first, Position(-99, -99, -99)); + fVertexStdevMap->emplace(clusterpair.first, -99); + } + }// end loop over the clusters +} + +//------------------------------------------------------------------------------ +std::vector VertexLeastSquares::FilterHits(std::vector hits) +{ + + // we can do this in two ways: + // 1. Filter the hits within a window to obtain a set of the most compatible hits with the first hit (fCompatibleHits = true) + // 2. Filter hits by checking if secondary hits are within 10ns of the mean of the first 4 hits (fCompatibleHits = false) + + // either way the idea is to shave off hits that do not come from the *assumed* point interaction / cherenkov cone (remove hits that could have come from reflections / secondary interactions) + + // if enabled (fGoodTimes == True) we can also reject hits that have poor timing std from the laser analysis + + std::vector filtered; + std::vector reliableHits; + + + // first, (if applicable) apply a filter to remove unreliable PMTs + if (fGoodTimes) { + for (auto& h : hits) { + int pmtid = h.GetTubeId(); + + // look up timing std in map + double timing_sigma = 1.0; // default + auto it = ChannelKeyToTimingSigmaMap->find(pmtid); + if (it != ChannelKeyToTimingSigmaMap->end()) { + timing_sigma = it->second; + logmessage = "VertexLeastSquares: PMTID " + std::to_string(pmtid) + ", timing std = " + std::to_string(timing_sigma) + "ns"; + Log(logmessage, v_debug, verbosity); + } else { + logmessage = "VertexLeastSquares: Didn't find timing std for PMTID " + std::to_string(pmtid) + ", using default = " + std::to_string(timing_sigma) + "ns"; + Log(logmessage, v_warning, verbosity); + } + + // reject unreliable PMTs + if (timing_sigma > fMinGoodTime) { + logmessage = "VertexLeastSquares: Rejecting PMT " + std::to_string(pmtid) + ", with timing std = " + std::to_string(timing_sigma) + "ns (threshold " + std::to_string(fMinGoodTime) + "ns)"; + Log(logmessage, v_debug, verbosity); + continue; + } + reliableHits.push_back(h); // keep hit if PMT std is low enough + } + + // bail if not enough reliable hits + if (reliableHits.size() < 4) { + Log("VertexLeastSquares: Not enough reliable hits (" + + std::to_string(reliableHits.size()) + ")", + v_warning, verbosity); + return filtered; // empty, RunLoop function will skip the reconstruction and assign -99 for the vertex + } + } + + + // pick which hits to use for the main filtering step + std::vector& hitsToFilter = fGoodTimes ? reliableHits : hits; + + + if (fCompatibleHits) { // option 1 + + // Sort the hits by time + std::sort(hitsToFilter.begin(), hitsToFilter.end(), + [](const Hit &a, const Hit &b) { return a.GetTime() < b.GetTime(); }); + + // Check for causally independent hits + // time difference between the hits must be less than the relative time of flight between the PMTs + std::map> compatibleIds; + int mostCompatibleIds = 0; + int bestIdx = 0; + for (uint idx = 0; idx < hitsToFilter.size(); ++idx) { + Detector *det_i = fGeom->ChannelToDetector(hitsToFilter[idx].GetTubeId()); + Position pos_i = det_i->GetDetectorPosition(); + + for (uint jdx = 0; jdx < hitsToFilter.size(); ++jdx) { + Detector *det_j = fGeom->ChannelToDetector(hitsToFilter[jdx].GetTubeId()); + Position pos_j = det_j->GetDetectorPosition(); + + double deltaT = abs(hitsToFilter[jdx].GetTime() - hitsToFilter[idx].GetTime()); + double tof = (pos_i - pos_j).Mag() / fSoL; + + if (deltaT <= tof) + compatibleIds[idx].insert(jdx); + } + + if (compatibleIds[idx].size() > mostCompatibleIds) { + mostCompatibleIds = compatibleIds[idx].size(); + bestIdx = idx; + } + } + + // Keep the set of most compatible ids if they're within 13 ns (light travel time across the inner structure) of the first hit + int firstHitIdx = *compatibleIds[bestIdx].begin(); + double t0 = hitsToFilter[firstHitIdx].GetTime(); + for (auto idx : compatibleIds[bestIdx]) { + if (abs(hitsToFilter[idx].GetTime() - t0) > 13) continue; + filtered.push_back(hitsToFilter[idx]); + } + + } else { // option 2 + + // Sort the hits by time + std::sort(hitsToFilter.begin(), hitsToFilter.end(), + [](const Hit &a, const Hit &b) { return a.GetTime() < b.GetTime(); }); + + // Keep the first four and find the mean time + double mean = 0; + for (uint idx = 0; idx < 4; ++idx) { + mean += hits[idx].GetTime(); + filtered.push_back(hitsToFilter[idx]); + } + mean = mean/4.; + + // Look through the rest of the hits, + // check that they are within 10 ns (for reflections), + // and check that they could come from the same point + for (uint idx = 4; idx < hitsToFilter.size(); ++idx) { + if ((hitsToFilter[idx].GetTime() - mean) > 10) + continue; + + bool keep = true; + for (uint jdx = 0; jdx < 4; ++jdx) { + Detector *det_i = fGeom->ChannelToDetector(hitsToFilter[idx].GetTubeId()); + Position pos_i = det_i->GetDetectorPosition(); + + Detector *det_j = fGeom->ChannelToDetector(filtered[jdx].GetTubeId()); + Position pos_j = det_j->GetDetectorPosition(); + + double deltaT = abs(filtered[jdx].GetTime() - hitsToFilter[idx].GetTime()); + double tof = (pos_i - pos_j).Mag() / fSoL; + + if (deltaT > tof) { + keep = false; + break; + } + } + + if (!keep) continue; + filtered.push_back(hitsToFilter[idx]); + } + + } + + return filtered; + +} + + +//------------------------------------------------------------------------------ +std::vector VertexLeastSquares::FilterHitsMC(std::vector hits) +{ + + std::vector filtered; + std::vector reliableHits; + + // (if applicable) apply a filter to remove unreliable PMTs + if (fGoodTimes) { + for (auto& h : hits) { + int pmtid = h.GetTubeId(); + + // look up timing std in map + double timing_sigma = 1.0; // default + auto it = ChannelKeyToTimingSigmaMap->find(pmtid); + if (it != ChannelKeyToTimingSigmaMap->end()) { + timing_sigma = it->second; + logmessage = "VertexLeastSquares: PMTID " + std::to_string(pmtid) + ", timing std = " + std::to_string(timing_sigma) + "ns"; + Log(logmessage, v_debug, verbosity); + } else { + logmessage = "VertexLeastSquares: Didn't find timing std for PMTID " + std::to_string(pmtid) + ", using default = " + std::to_string(timing_sigma) + "ns"; + Log(logmessage, v_warning, verbosity); + } + + // reject unreliable PMTs + if (timing_sigma > fMinGoodTime) { + logmessage = "VertexLeastSquares: Rejecting PMT " + std::to_string(pmtid) + ", with timing std = " + std::to_string(timing_sigma) + "ns (threshold " + std::to_string(fMinGoodTime) + "ns)"; + Log(logmessage, v_debug, verbosity); + continue; + } + reliableHits.push_back(h); // keep hit if PMT std is low enough + } + + // bail if not enough reliable hits + if (reliableHits.size() < 4) { + Log("VertexLeastSquares: Not enough reliable hits (" + + std::to_string(reliableHits.size()) + ")", + v_warning, verbosity); + return filtered; // empty, RunLoop function will skip the reconstruction and assign -99 for the vertex + } + } + + // pick which hits to use for the main filtering step + std::vector& hitsToFilter = fGoodTimes ? reliableHits : hits; + + if (fCompatibleHits) { // option 1 + + // Sort the hits by time + std::sort(hitsToFilter.begin(), hitsToFilter.end(), + [](const MCHit &a, const MCHit &b) {return a.GetTime() < b.GetTime();}); + + // Check for causally independent hits + // time difference between the hits must be less than the relative time of flight between the PMTs + std::map> compatibleIds; + int mostCompatibleIds = 0; + int bestIdx = 0; + for (uint idx = 0; idx < hitsToFilter.size(); ++idx) { + Detector *det_i = fGeom->ChannelToDetector(hitsToFilter[idx].GetTubeId()); + Position pos_i = det_i->GetDetectorPosition(); + + for (uint jdx = 0; jdx < hitsToFilter.size(); ++jdx) { + Detector *det_j = fGeom->ChannelToDetector(hitsToFilter[jdx].GetTubeId()); + Position pos_j = det_j->GetDetectorPosition(); + + double deltaT = abs(hitsToFilter[jdx].GetTime() - hitsToFilter[idx].GetTime()); + double tof = (pos_i - pos_j).Mag() / fSoL; + + if (deltaT <= tof) + compatibleIds[idx].insert(jdx); + } + + if (compatibleIds[idx].size() > mostCompatibleIds) { + mostCompatibleIds = compatibleIds[idx].size(); + bestIdx = idx; + } + } + + // Keep the set of most compatible ids if they're within 13 ns (light travel time across the inner structure) of the first hit + int firstHitIdx = *compatibleIds[bestIdx].begin(); + double t0 = hitsToFilter[firstHitIdx].GetTime(); + for (auto idx : compatibleIds[bestIdx]) { + if (abs(hitsToFilter[idx].GetTime() - t0) > 13) continue; + filtered.push_back(hitsToFilter[idx]); + } + + } else { + + // Sort the hits by time + std::sort(hitsToFilter.begin(), hitsToFilter.end(), + [](const MCHit &a, const MCHit &b) {return a.GetTime() < b.GetTime();}); + + // Keep the first four and find the mean time + double mean = 0; + for (uint idx = 0; idx < 4; ++idx) { + mean += hitsToFilter[idx].GetTime(); + filtered.push_back(hitsToFilter[idx]); + } + mean = mean/4.; + + // Look through the rest of the hits, + // check that they are within 10 ns (for reflections), + // and check that they could come from the same point + for (uint idx = 4; idx < hitsToFilter.size(); ++idx) { + if ((hitsToFilter[idx].GetTime() - mean) > 10) + continue; + + bool keep = true; + for (uint jdx = 0; jdx < 4; ++jdx) { + Detector *det_i = fGeom->ChannelToDetector(hitsToFilter[idx].GetTubeId()); + Position pos_i = det_i->GetDetectorPosition(); + + Detector *det_j = fGeom->ChannelToDetector(filtered[jdx].GetTubeId()); + Position pos_j = det_j->GetDetectorPosition(); + + double deltaT = abs(filtered[jdx].GetTime() - hitsToFilter[idx].GetTime()); + double tof = (pos_i - pos_j).Mag() / fSoL; + + if (deltaT > tof) { + keep = false; + break; + } + } + + if (!keep) continue; + filtered.push_back(hitsToFilter[idx]); + } + + } + + return filtered; +} diff --git a/UserTools/VertexLeastSquares/VertexLeastSquares.h b/UserTools/VertexLeastSquares/VertexLeastSquares.h new file mode 100644 index 000000000..a399d9438 --- /dev/null +++ b/UserTools/VertexLeastSquares/VertexLeastSquares.h @@ -0,0 +1,109 @@ +#ifndef VertexLeastSquares_H +#define VertexLeastSquares_H + +#include +#include + +#include "Tool.h" +#include "MatrixUtil.h" +#include "Position.h" +#include "Hit.h" + +#include "TFile.h" +#include "TTree.h" + + + +/** + * \class VertexLeastSquares + * +* +* $Author(s): S.Doran + A.Sutton$ +* $Date: 2025/09/26 00:00:00 $ +* Contact: doran@iastate.edu +*/ +class VertexLeastSquares: public Tool { + + + public: + + VertexLeastSquares(); ///< Simple constructor + bool Initialise(std::string configfile,DataModel &data); ///< Initialise Function for setting up Tool resources. @param configfile The path and name of the dynamic configuration file to read in. @param data A reference to the transient data class used to pass information between Tools. + bool Execute(); ///< Execute function used to perform Tool purpose. + bool Finalise(); ///< Finalise function used to clean up resources. + + // Generate a cylinder with XYZ nodes in each dimension + std::vector GenerateVetices(); + + // Evaluate the Jacobian and the solution vector at the guess vertex + // we'll be solving A*dx = b where A = Jacobian(guess), b = -f(guess), + // f = vector of the hit time functions: 0 = 1/c * sqrt( (X-xi)^2 + (Y-yi)^2 + (Z-zi)^2 ) + T - ti + // X,Y,Z are the coords of the vertex and T is the emission time of the light + // xi,yi,zi are the PMT coords and ti is the PMT hit time + void EvalAtGuessVertex(util::Matrix &A, util::Vector &b, const Position &guess, const std::vector &hits); + void EvalAtGuessVertexMC(util::Matrix &A, util::Vector &b, const Position &guess, const std::vector &hits); + + // Actually run the thing until convergence + void RunLoop(); + void RunLoopMC(); + + // Filter hits based on time separation and causality + std::vector FilterHits(std::vector hits); + std::vector FilterHitsMC(std::vector hits); + + + private: + + // Configuration parameters + bool fUseMCHits; + double fBreakDist; + bool fDebugTree; + double fYSpacing; + int fNPlanarPoints; + double fRegularizer; + int fFittingSteps; + bool fExternalSeeding; + double fYBuffer; + double fRadialBuffer; + double fExternalVertexMaxRadius; + double fExternalVertexMaxHeight; + bool fCompatibleHits; + bool fGoodTimes; + double fMinGoodTime; + + // The ANNIE geometry service + Geometry *fGeom = nullptr; + + // timing uncertainty map + std::map* ChannelKeyToTimingSigmaMap; + + // The clusters we'll load from the CStore + std::map> *fClusterMap = nullptr; + std::map> *fClusterMapMC = nullptr; + + // Vertices and stdev we'll save to the ANNIEEvent + std::map *fVertexMap = nullptr; + std::map *fVertexStdevMap = nullptr; + + // The number of boundary points for the XZ plane verticies + int fNBoundary; + + // Golden ratio and pi + const double phi = (1. + sqrt(5.))/2.; + const double phisq = pow(phi, 2); + const double pi = 3.141592; + + // Speed of light in water, m/ns + const double fSoL = 0.299792458 * 3/4; + + /// \brief verbosity levels: if 'verbosity' < this level, the message type will be logged. + int verbosity; + int v_error=0; + int v_warning=1; + int v_message=2; + int v_debug=3; + std::string logmessage; +}; + + +#endif