[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFmHelix.h

///////////////////////////////////////////////////////////////////////////
//                                                                       //
// AliFmHelix: a helper helix class                                      //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#ifndef ALIFMHELIX_H
#define ALIFMHELIX_H

#include <math.h>
#include <utility>
#include <algorithm>
#include "AliFmThreeVector.h"
#if !defined(ST_NO_NAMESPACES)
using std::pair;
using std::swap;
using std::max;
#endif

#ifdef WIN32
#include "gcc2vs.h"
#endif

#ifdef __SUNPRO_CC
#include <ieeefp.h>
#define __FUNCTION__ "__FILE__:__LINE__"
#endif

class AliFmHelix {
public:
    /// curvature, dip angle, phase, origin, h
    AliFmHelix(double c, double dip, double phase,
	       const AliFmThreeVector<double>& o, int h=-1);
    
    virtual ~AliFmHelix();
    // AliFmHelix(const AliFmHelix&);			// use default
    // AliFmHelix& operator=(const AliFmHelix&);	// use default

    double       DipAngle()   const;           
    double       Curvature()  const;	/// 1/R in xy-plane
    double       Phase()      const;	/// aziumth in xy-plane measured from ring center
    double       XCenter()    const;	/// x-center of circle in xy-plane
    double       YCenter()    const;	/// y-center of circle in xy-plane
    int          H()          const;	/// -sign(q*B);
    
    const AliFmThreeVector<double>& Origin() const;	/// starting point

    void SetParameters(double c, double dip, double phase, const AliFmThreeVector<double>& o, int h);
    
    double       X(double s)  const;
    double       Y(double s)  const;
    double       Z(double s)  const;

    AliFmThreeVector<double>  At(double s) const;

    /// returns period length of helix
    double       Period()       const;
    
    /// path length at given r (cylindrical r)
    pair<double, double> PathLength(double r)   const;
    
    /// path length at given r (cylindrical r, cylinder axis at x,y)
    pair<double, double> PathLength(double r, double x, double y, bool scanPeriods = true);
    
    /// path length at distance of closest approach to a given point
    double       PathLength(const AliFmThreeVector<double>& p, bool scanPeriods = true) const;
    
    /// path length at intersection with plane
    double       PathLength(const AliFmThreeVector<double>& r,
			    const AliFmThreeVector<double>& n) const;

    /// path length at distance of closest approach in the xy-plane to a given point
    double       PathLength(double x, double y) const;

    /// path lengths at dca between two helices 
    pair<double, double> PathLengths(const AliFmHelix& h, bool scanPeriods = true) const;
    
    /// minimal distance between point and helix
    double       Distance(const AliFmThreeVector<double>& p, bool scanPeriods = true) const;    
    
    /// checks for valid parametrization
    bool         Valid(double world = 1.e+5) const {return !Bad(world);}
    int            Bad(double world = 1.e+5) const;
    
    /// move the origin along the helix to s which becomes then s=0
    virtual void MoveOrigin(double s);
    
    static const double fgkNoSolution;  // coinstant indicating lack of solution
    
protected:
    AliFmHelix();
    
    void SetCurvature(double d);	/// performs also various checks   
    void SetPhase(double d);	        
    void SetDipAngle(double d);
    
    /// value of S where distance in x-y plane is minimal
    double FudgePathLength(const AliFmThreeVector<double>& v) const;
    
protected:
    bool                   fSingularity;	// true for straight line case (B=0)
    AliFmThreeVector<double>  fOrigin;          // Helix origin
    double                 fDipAngle;           // Helix dip angle
    double                 fCurvature;          // curvature
    double                 fPhase;              // phase
    int                    fH;			// -sign(q*B);

    double                 fCosDipAngle;        // cosine of the dip angle
    double                 fSinDipAngle;        // sine of the dip angle
    double                 fCosPhase;           // cosine of the phase
    double                 fSinPhase;           // sine of the phase
#ifdef __ROOT__
  ClassDef(AliFmHelix,1)
#endif
};

//
//     Non-member functions
//
int operator== (const AliFmHelix&, const AliFmHelix&);
int operator!= (const AliFmHelix&, const AliFmHelix&);
ostream& operator<<(ostream&, const AliFmHelix&);

//
//     Inline functions
//
inline int AliFmHelix::H() const {return fH;}

inline double AliFmHelix::DipAngle() const {return fDipAngle;}

inline double AliFmHelix::Curvature() const {return fCurvature;}

inline double AliFmHelix::Phase() const {return fPhase;}

inline double AliFmHelix::X(double s) const
{
    if (fSingularity)
	return fOrigin.x() - s*fCosDipAngle*fSinPhase;
    else
	return fOrigin.x() + (cos(fPhase + s*fH*fCurvature*fCosDipAngle)-fCosPhase)/fCurvature;
}
 
inline double AliFmHelix::Y(double s) const
{
    if (fSingularity)
	return fOrigin.y() + s*fCosDipAngle*fCosPhase;
    else
	return fOrigin.y() + (sin(fPhase + s*fH*fCurvature*fCosDipAngle)-fSinPhase)/fCurvature;
}

inline double AliFmHelix::Z(double s) const
{
    return fOrigin.z() + s*fSinDipAngle;
}

inline const AliFmThreeVector<double>& AliFmHelix::Origin() const {return fOrigin;}

inline AliFmThreeVector<double> AliFmHelix::At(double s) const
{
    return AliFmThreeVector<double>(X(s), Y(s), Z(s));
}

inline double AliFmHelix::PathLength(double x, double y) const
{
    return FudgePathLength(AliFmThreeVector<double>(x, y, 0));
}
inline int AliFmHelix::Bad(double WorldSize) const
{

    int ierr;
    if (!::finite(fDipAngle    )) 	return   11;
    if (!::finite(fCurvature   )) 	return   12;

    ierr = fOrigin.Bad(WorldSize);
    if (ierr)                           return    3+ierr*100;

    if (::fabs(fDipAngle)  >1.58)	return   21;
    double qwe = ::fabs(::fabs(fDipAngle)-M_PI/2);
    if (qwe < 1./WorldSize      ) 	return   31; 

    if (::fabs(fCurvature) > WorldSize)	return   22;
    if (fCurvature < 0          )	return   32;

    if (abs(fH) != 1            )       return   24; 

    return 0;
}

#endif
Commit	Line	Data
d0e92d9a	1	///////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFmHelix: a helper helix class //
	4	// //
	5	///////////////////////////////////////////////////////////////////////////
	6
	7	#ifndef ALIFMHELIX_H
	8	#define ALIFMHELIX_H
	9
	10	#include <math.h>
	11	#include <utility>
	12	#include <algorithm>
	13	#include "AliFmThreeVector.h"
	14	#if !defined(ST_NO_NAMESPACES)
	15	using std::pair;
	16	using std::swap;
	17	using std::max;
	18	#endif
	19
	20	#ifdef WIN32
	21	#include "gcc2vs.h"
	22	#endif
	23
6ac81678	24	#ifdef __SUNPRO_CC
	25	#include <ieeefp.h>
	26	#define __FUNCTION__ "__FILE__:__LINE__"
	27	#endif
	28
d0e92d9a	29	class AliFmHelix {
	30	public:
	31	/// curvature, dip angle, phase, origin, h
	32	AliFmHelix(double c, double dip, double phase,
	33	const AliFmThreeVector<double>& o, int h=-1);
	34
	35	virtual ~AliFmHelix();
	36	// AliFmHelix(const AliFmHelix&); // use default
	37	// AliFmHelix& operator=(const AliFmHelix&); // use default
	38
	39	double DipAngle() const;
	40	double Curvature() const; /// 1/R in xy-plane
	41	double Phase() const; /// aziumth in xy-plane measured from ring center
	42	double XCenter() const; /// x-center of circle in xy-plane
	43	double YCenter() const; /// y-center of circle in xy-plane
	44	int H() const; /// -sign(q*B);
	45
	46	const AliFmThreeVector<double>& Origin() const; /// starting point
	47
	48	void SetParameters(double c, double dip, double phase, const AliFmThreeVector<double>& o, int h);
	49
	50	double X(double s) const;
	51	double Y(double s) const;
	52	double Z(double s) const;
	53
	54	AliFmThreeVector<double> At(double s) const;
	55
	56	/// returns period length of helix
	57	double Period() const;
	58
	59	/// path length at given r (cylindrical r)
	60	pair<double, double> PathLength(double r) const;
	61
	62	/// path length at given r (cylindrical r, cylinder axis at x,y)
	63	pair<double, double> PathLength(double r, double x, double y, bool scanPeriods = true);
	64
	65	/// path length at distance of closest approach to a given point
	66	double PathLength(const AliFmThreeVector<double>& p, bool scanPeriods = true) const;
	67
	68	/// path length at intersection with plane
	69	double PathLength(const AliFmThreeVector<double>& r,
	70	const AliFmThreeVector<double>& n) const;
	71
	72	/// path length at distance of closest approach in the xy-plane to a given point
	73	double PathLength(double x, double y) const;
	74
	75	/// path lengths at dca between two helices
	76	pair<double, double> PathLengths(const AliFmHelix& h, bool scanPeriods = true) const;
	77
	78	/// minimal distance between point and helix
	79	double Distance(const AliFmThreeVector<double>& p, bool scanPeriods = true) const;
	80
	81	/// checks for valid parametrization
	82	bool Valid(double world = 1.e+5) const {return !Bad(world);}
	83	int Bad(double world = 1.e+5) const;
	84
	85	/// move the origin along the helix to s which becomes then s=0
	86	virtual void MoveOrigin(double s);
	87
d92ed900	88	static const double fgkNoSolution; // coinstant indicating lack of solution
d0e92d9a	89
	90	protected:
	91	AliFmHelix();
	92
	93	void SetCurvature(double d); /// performs also various checks
	94	void SetPhase(double d);
	95	void SetDipAngle(double d);
	96
	97	/// value of S where distance in x-y plane is minimal
	98	double FudgePathLength(const AliFmThreeVector<double>& v) const;
	99
	100	protected:
	101	bool fSingularity; // true for straight line case (B=0)
	102	AliFmThreeVector<double> fOrigin; // Helix origin
	103	double fDipAngle; // Helix dip angle
	104	double fCurvature; // curvature
	105	double fPhase; // phase
	106	int fH; // -sign(q*B);
	107
	108	double fCosDipAngle; // cosine of the dip angle
	109	double fSinDipAngle; // sine of the dip angle
	110	double fCosPhase; // cosine of the phase
	111	double fSinPhase; // sine of the phase
	112	#ifdef __ROOT__
	113	ClassDef(AliFmHelix,1)
	114	#endif
	115	};
	116
	117	//
	118	// Non-member functions
	119	//
	120	int operator== (const AliFmHelix&, const AliFmHelix&);
	121	int operator!= (const AliFmHelix&, const AliFmHelix&);
	122	ostream& operator<<(ostream&, const AliFmHelix&);
	123
	124	//
	125	// Inline functions
	126	//
	127	inline int AliFmHelix::H() const {return fH;}
	128
	129	inline double AliFmHelix::DipAngle() const {return fDipAngle;}
	130
	131	inline double AliFmHelix::Curvature() const {return fCurvature;}
	132
	133	inline double AliFmHelix::Phase() const {return fPhase;}
	134
	135	inline double AliFmHelix::X(double s) const
	136	{
	137	if (fSingularity)
	138	return fOrigin.x() - sfCosDipAnglefSinPhase;
	139	else
	140	return fOrigin.x() + (cos(fPhase + sfHfCurvature*fCosDipAngle)-fCosPhase)/fCurvature;
	141	}
	142
	143	inline double AliFmHelix::Y(double s) const
	144	{
	145	if (fSingularity)
	146	return fOrigin.y() + sfCosDipAnglefCosPhase;
	147	else
	148	return fOrigin.y() + (sin(fPhase + sfHfCurvature*fCosDipAngle)-fSinPhase)/fCurvature;
	149	}
	150
	151	inline double AliFmHelix::Z(double s) const
	152	{
153	return fOrigin.z() + s*fSinDipAngle;
154	}
155
156	inline const AliFmThreeVector<double>& AliFmHelix::Origin() const {return fOrigin;}
157
158	inline AliFmThreeVector<double> AliFmHelix::At(double s) const
159	{
160	return AliFmThreeVector<double>(X(s), Y(s), Z(s));
161	}
162
163	inline double AliFmHelix::PathLength(double x, double y) const
164	{
165	return FudgePathLength(AliFmThreeVector<double>(x, y, 0));
166	}
167	inline int AliFmHelix::Bad(double WorldSize) const
168	{
169
170	int ierr;
171	if (!::finite(fDipAngle )) return 11;
172	if (!::finite(fCurvature )) return 12;
173
a19edcc9	174	ierr = fOrigin.Bad(WorldSize);
d0e92d9a	175	if (ierr) return 3+ierr*100;
	176
	177	if (::fabs(fDipAngle) >1.58) return 21;
	178	double qwe = ::fabs(::fabs(fDipAngle)-M_PI/2);
	179	if (qwe < 1./WorldSize ) return 31;
	180
	181	if (::fabs(fCurvature) > WorldSize) return 22;
	182	if (fCurvature < 0 ) return 32;
	183
	184	if (abs(fH) != 1 ) return 24;
	185
	186	return 0;
	187	}
	188
	189	#endif