Template specialisation of fit function of acoustic model based on linear approximation. More...

#include <JKatoomba.hh>

Inheritance diagram for JACOUSTICS::JKatoomba< JEstimator >:

Public Member Functions
	JKatoomba (const JGeometry &geometry, const JSoundVelocity &velocity, const int option)
	Constructor. More...

template<class T >
const JModel &	operator() (T __begin, T __end) const
	Get start values of string parameters. More...

Static Public Attributes
static constexpr double	TOLERANCE = 1.0e-8
	Tolerance for SVD. More...

Private Attributes
JModel	value

Detailed Description

template<>
struct JACOUSTICS::JKatoomba< JEstimator >

Template specialisation of fit function of acoustic model based on linear approximation.

Definition at line 334 of file JKatoomba.hh.

Constructor & Destructor Documentation

JACOUSTICS::JKatoomba< JEstimator >::JKatoomba	(	const JGeometry &	geometry,
		const JSoundVelocity &	velocity,
		const int	option
	)

inline

Constructor.

Parameters

geometry	detector geometry
velocity	sound velocity
option	option

Definition at line 346 of file JKatoomba.hh.

                                             :
       JKatoomba<>(geometry, velocity, option)
     {};

Member Function Documentation

template<class T >

const JModel& JACOUSTICS::JKatoomba< JEstimator >::operator()	(	T	__begin,
		T	__end
	)		const

inline

Get start values of string parameters.

Note that this method may throw an exception.

Parameters

__begin	begin of hits
__end	end of hits

Returns: model

Definition at line 362 of file JKatoomba.hh.

     {
       using namespace std;
       using namespace JPP;
       using namespace JGEOMETRY;
 
 #ifdef THREAD_SAFE
       static mutex mtx;                               // TDecompSVD
 
       unique_lock<mutex> lock(mtx);
 #endif
       value = JModel(__begin, __end);                 // set up global model according data
 
       if (this->option == JMODEL::fit_emitters_only_t)
         value.setOption(JMODEL::fit_emitters_only_t);
       else
         value.setOption(JMODEL::fit_emitters_and_strings_1st_order_t);
 
       H_t H;                                          // H-equation as per hit
       I_t i;                                          // indices of H-equation in global model 
 
       const size_t N = value.getN();
 
       TMatrixD       V(N,N);                          // V = H^T x H
       vector<double> Y(N, 0.0);                       // Y = H^T x y
 
       for (T hit = __begin; hit != __end; ++hit) {
 
         const JString&    string   = geometry[hit->getString()];
         const JPosition3D position = string.getPosition(hit->getFloor());
         const double      Vi       = velocity.getInverseVelocity(hit->getDistance(position), hit->getZ(), position.getZ());
 
         const double      h1 = string.getHeight(hit->getFloor());
         const JPosition3D p1 = string.getPosition()  -  hit->getPosition();
         const double      ds = sqrt(p1.getLengthSquared()  +  h1*h1  +  2.0*p1.getZ()*h1);
         const double      y  = hit->getValue()  -  Vi*ds;
         const double      W  = sqrt(hit->getWeight());
 
         H.t1  =  W * 1.0;
         H.tx  =  W * Vi * p1.getX() * h1 / ds;
         H.ty  =  W * Vi * p1.getY() * h1 / ds;
         
         i.t1  =  value.getIndex(hit->getEKey(),   &H_t::t1);
         i.tx  =  value.getIndex(hit->getString(), &H_t::tx);
         i.ty  =  value.getIndex(hit->getString(), &H_t::ty);
 
         V(i.t1, i.t1) += H.t1 * H.t1;
 
         Y[i.t1] += W * H.t1 * y;
 
         if (hit->getFloor() != 0) {
 
           if (this->option != JMODEL::fit_emitters_only_t) {
           
             V(i.t1, i.tx) += H.t1 * H.tx;  V(i.t1, i.ty) += H.t1 * H.ty;
             V(i.tx, i.t1) += H.tx * H.t1;  V(i.ty, i.t1) += H.ty * H.t1;
 
             V(i.tx, i.tx) += H.tx * H.tx;  V(i.tx, i.ty) += H.tx * H.ty;
             V(i.ty, i.tx) += H.ty * H.tx;  V(i.ty, i.ty) += H.ty * H.ty;
 
             Y[i.tx] += W * H.tx * y;
             Y[i.ty] += W * H.ty * y;
           }
         }
       }
 
       // evaluate (H^T x H)^-1 x H^T x y 
 
       TDecompSVD svd(V, TOLERANCE);
 
       Bool_t status;
 
       V = svd.Invert(status);
 
       for (size_t row = 0; row != N; ++row) {
         for (size_t col = 0; col != N; ++col) {
           value[row] += V(row,col) * Y[col];
         }
       }
       
       return value;
     }