ElasticArms.h

Go to the documentation of this file.

#ifndef EARMS_ELASTICARMS_H
#define EARMS_ELASTICARMS_H
#include <vector>

namespace earms {
  ///
  /// Fit a vertex and prongs to a collection of hits using the
  /// "elastic arms" method
  ///
  class ElasticArms {
  public:
    ///
    /// Construct an elastic arms fitter
    ///
    /// \param n       - number of hits
    /// \param m       - number of prongs
    /// \param lambda  - noise penalty
    /// \param lambdav - distance scale from vertex to first hit
    ///
    ElasticArms(unsigned int n, unsigned int m, double lambda, double lambdav);
    
    ///
    /// Add a hit to the fit
    ///
    /// \param i     - hit index [0...N)
    /// \param z     - z location of hit (cm)
    /// \param xory  - x or y location of hit (cm)
    /// \param sigma - hit location uncertainty (cm)
    /// \param view  - X view or Y view?
    ///
    void SetHit(unsigned int i, double z, double xory, double sigma, int view);

    ///
    /// The the annealing temperature.
    ///
    /// \param t - the temperature t is in units of chi^2
    ///
    void SetT(double t);

    ///
    /// The the noise penalty
    ///
    /// \param lam - the noise penalty in units of chi^2
    ///
    void SetLambda(double lam);
    
    ///
    /// Update the matrix of associations
    ///
    void UpdateVia();

    ///
    /// Update the vector of arm-vertex distances
    ///
    void UpdateDaa();
    
    ///
    /// Set the vertex parameters
    ///
    /// \param  x - x position in cm
    /// \param  y - y position in cm
    /// \param  z - z position in cm
    ///
    void SetVertex(double x, double y, double z);

    void GetVertex(double& x, double& y, double& z) const;

    ///
    /// Set the arm direction parameters
    ///
    /// \param theta - zenith angle in radians
    /// \param phi   - azimuthal angle in radians
    ///
    void SetArm(unsigned int i, double theta, double phi);
    
    ///
    /// Compute the distances^2 from hit i to track a
    ///
    void UpdateMia();
    
    ///
    /// Return the "energy" function, a chi^2-like goodness of fit
    ///
    double E();

    ///
    /// Calculation the total weight of an arm a. Roughly speaking
    /// this is the number of hits on the arm.
    ///
    double SumVi(unsigned int a) const;

    ///
    /// Find the leading arm. The with the highest weight
    ///
    unsigned int LeadingArm() const;
    
  private:
    int    fProngFit;      ///< 1 - treat arms as prongs in fit, 0 - treat as lines
    double fT;             ///< A measure of the "temperature" in units of sigma^2
    double fBeta;          ///< 1/T in units of 1/sigma^2
    bool   fBetaUnchanged; ///< Whether fBeta is the same since last call to UpdateMia
    double fLambda;        ///< Noise penalty term in units of sigma^2
    double fExpBetaLambda; ///< exp(-fBeta*fLambda) - computed once used many
    double fLambdaV;       ///< Distance scale for hits to appear on (cm)

  public:
    /// Data to fit
    unsigned int        fN;        ///< Number of hits
    std::vector<double> fXorY;     ///< Transverse hit positions (cm)
    std::vector<int>    fView;     ///< Which view is the hit in?
    std::vector<double> fZ;        ///< Z locations of hits (cm)
    std::vector<double> fWsqr;     ///< Combined weight for hit 1/(sigmaH^2 + sigmaMS^2)
    std::vector<double> fDhi;      ///< From hit to vertex location^2 [cm^2]
    std::vector<double> fDaa;      ///< Distance from vertex to prong [cm]

    typedef std::vector< std::vector<double> > matrix;
    matrix fMia;        ///< Distance^2 hit i to track a (sigma^2)
    std::vector< std::vector<bool> > fMiaUnchanged; ///< Is fMia changed in UpdateMia?
    std::vector<bool> fMiaRowUnchanged; ///< Is a row of fMia changed in UpdateMia?
    std::vector<bool> fMiaColUnchanged; ///< Is a column of fMia changed in UpdateMia?
    matrix fExpBetaMia; ///< exp(-fBeta*Mia) - compute once used many time
    matrix fVia;        ///< Strength of association hit i to track a
    
  public:
    /// Model parameters
    double              fX0;       ///< Vertex x location (cm)
    double              fY0;       ///< Vertex y location (cm)
    double              fZ0;       ///< Vertex z location (cm)
    unsigned int        fM;        ///< Number of arms
    std::vector<double> fdXds; ///< Track dx/ds
    std::vector<double> fdYds; ///< Track dy/ds
    std::vector<double> fdZds; ///< Track dz/ds
  };
}

#endif
////////////////////////////////////////////////////////////////////////