include/IRT2/IRT.h

0001 #pragma once
0002
0003 #include <vector>
0004
0005 #include <TVector3.h>
0006 #include <TObject.h>
0007 #include <TRef.h>
0008
0009 // Don't even care about the units here; [mm] in fact; just make it practically small;
0010 #define _IRT_PRECISION_DEFAULT_     (1E-7)
0011 // Well, assume radians in [theta,phi] -> 1urad step should be good enough;
0012 #define _IRT_JACOBIAN_STEP_DEFAULT_ (1E-6)
0013 #define _IRT_ITERATION_LIMIT_        ( 50)
0014 // Here perhaps want something tangible, on a [mm] scale (?), matching the uncertaintly
0015 // in the sensor position measurement, and emission point uncertainty along the trajectory;
0016 #define _IRT_DERIVATIVE_XYZ_STEP_    (1.00)
0017 // Refractive index variation; in units of the (1-n); ~1% makes sense?;
0018 //#define _IRT_DERIVATIVE_NNN_STEP_    (0.01)
0019
0020 #include "CherenkovRadiator.h"
0021 #include "OpticalBoundary.h"
0022 #include "ParametricSurface.h"
0023 #include "IRTSolution.h"
0024
0025 namespace IRT2 {
0026
0027 class IRT: public TObject {
0028  public:
0029   IRT(/*unsigned sector = 0*/): /*m_Sector(sector),*/ m_IterationLimit(_IRT_ITERATION_LIMIT_),
0030     m_Precision(_IRT_PRECISION_DEFAULT_), m_JacobianStep(_IRT_JACOBIAN_STEP_DEFAULT_) {};
0031   ~IRT() {};
0032
0033   void AddOpticalBoundary(OpticalBoundary *boundary) {
0034     _m_OpticalBoundaries.push_back(boundary);
0035   };
0036
0037   void SetIterationLimit(unsigned value) { m_IterationLimit = value; };
0038   //void SetPrecision   (double value) { m_Precision = value; };
0039   void SetJacobianStep(double value) { m_JacobianStep = value; };
0040
0041   IRTSolution Solve(const TVector3 &xfrom, const TVector3 &nfrom, const TVector3 &xto, const TVector3 &beam,
0042             bool derivatives = false, const IRTSolution *seed = 0);
0043
0044   //unsigned GetSector( void ) const { return m_Sector; };
0045
0046   inline OpticalBoundary *tail( void ) const {
0047     return _m_OpticalBoundaries.size() ? GetOpticalBoundary(_m_OpticalBoundaries.size()-1) : 0;
0048   };
0049
0050  private:
0051   bool Transport(const TVector3 &xfrom, const TVector3 &nfrom, double *length = 0);
0052   inline OpticalBoundary *GetOpticalBoundary(unsigned id) const {
0053     return (id < _m_OpticalBoundaries.size() ?
0054             dynamic_cast<OpticalBoundary*>(_m_OpticalBoundaries[id].GetObject()) : 0);
0055   };
0056
0057   IRTSolution Solve(const TVector3 &xfrom, const TVector3 &nfrom, const double xy[2], const TVector3 &beam,
0058             bool derivatives = false, const IRTSolution *seed = 0);
0059
0060   // FIXME: this is not the right place for this variable;
0061   //unsigned m_Sector;
0062
0063   unsigned m_IterationLimit;
0064   double m_Precision, m_JacobianStep;
0065
0066   std::vector<TRef> _m_OpticalBoundaries;
0067
0068 #ifndef DISABLE_ROOT_IO
0069   ClassDef(IRT, 4);
0070 #endif
0071 };
0072
0073 } // namespace IRT2