MpdTofGeoUtils.h 7.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177
  1. //------------------------------------------------------------------------------------------------------------------------
  2. #ifndef __MPD_TOF_GEO_UTILS_H
  3. #define __MPD_TOF_GEO_UTILS_H 1
  4. #include <TVector3.h>
  5. #include <TList.h>
  6. #include <TGeoMatrix.h>
  7. #include "IntervalTree.h"
  8. //------------------------------------------------------------------------------------------------------------------------
  9. class LRectangle // convex quadrangle
  10. {
  11. bool IsInvalid;
  12. inline void getRPhi(const TVector3& point, Double_t& Rmin, Double_t& Rmax, Double_t& Phimin, Double_t& Phimax)
  13. {
  14. const Double_t R = point.Perp(), Phi = point.Phi();
  15. if(R < Rmin) Rmin = R;
  16. if(R > Rmax) Rmax = R;
  17. if(Phi < Phimin) Phimin = Phi;
  18. if(Phi > Phimax) Phimax = Phi;
  19. }
  20. public:
  21. enum Side_t { kRight=0, kLeft=1, kInvalid= -1 };
  22. Int_t volumeUID;
  23. TVector3 A, B, C, D, center, perp; // [cm]
  24. LRectangle() : IsInvalid(true), volumeUID(kInvalid) {};
  25. LRectangle(Int_t uid, const TVector3& a, const TVector3& b, const TVector3& c, const TVector3& d, bool check = false);
  26. TVector3 GetCenter() const{ return (A+B+C+D) * 0.25;}
  27. bool isInvalid() const{ return IsInvalid;}
  28. Double_t GetIntersectionSquare(const LRectangle& v)const;
  29. void GetRPhiRanges(Double_t& Rmin, Double_t& Rmax, Double_t& Phimin, Double_t& Phimax);
  30. Double_t DistanceFromPointToLineSegment(const TVector3* pos, const TVector3& P1,const TVector3& P2)const;
  31. Double_t DistanceFromPointToLine(const TVector3* pos, const TVector3& P1,const TVector3& P2)const;
  32. Double_t MinDistanceToEdge(const TVector3* pos, Side_t& side) const;
  33. Double_t Angle(const LRectangle& r)const{return perp.Angle(r.perp);}; // [rad]
  34. void Rotate(const TRotation& rot);
  35. void Shift(const TVector3& shift){ A +=shift; B +=shift; C +=shift; D +=shift; };
  36. void Dump(const char* comment = nullptr, std::ostream& out = std::cout) const;
  37. static void Test(const char* comment = nullptr, std::ostream& out = std::cout);
  38. //-----------------------------------------------------
  39. inline double Square()const
  40. {
  41. return (B - A).Cross(D - A).Mag();
  42. }
  43. //-----------------------------------------------------
  44. inline bool IsParallelPlane(const LRectangle& v)const
  45. {
  46. // return (perp == v.perp);
  47. const static double epsilon = 1.e-6;
  48. return ( (std::fabs(perp.X() - v.perp.X()) < epsilon) && (std::fabs(perp.Y() - v.perp.Y()) < epsilon) && (std::fabs(perp.Z() - v.perp.Z()) < epsilon) );
  49. }
  50. //-----------------------------------------------------
  51. inline std::pair<bool,double> IsSamePlane(const LRectangle& v)const
  52. {
  53. if(!IsParallelPlane(v)) return {false, 0.};
  54. const static double epsilon = 1.e-6;
  55. double D1 = A*perp, D2 = v.A* v.perp, distance = fabs(D2 - D1);
  56. if(distance < epsilon) return {true, 0.};
  57. return {false, distance};
  58. }
  59. //-----------------------------------------------------
  60. inline void InitCenterPerp()
  61. {
  62. center = (A+B+C+D) * 0.25;
  63. perp = (B-A).Cross(D-A).Unit();
  64. }
  65. //-----------------------------------------------------
  66. inline void CheckInValid()
  67. {
  68. IsInvalid = false;
  69. // Convex Polygon Definition: A polygon that has all interior angles less than 180°
  70. ;
  71. // Sum of Interior Angles, sum = 180*(n-2) degree, where n is the number of sides
  72. // A square has 4 sides, so interior angles sum = 360°
  73. ;
  74. // Rectangle check - all angles == 90 degree
  75. TVector3 ab = A-B, bc = B-C, cd = C-D, da = D-A;
  76. if( ab.Dot(bc) != 0. || bc.Dot(cd) != 0. || cd.Dot(da) != 0. || da.Dot(ab) != 0.)
  77. {
  78. std::cerr<<"\n ---> ERROR: invalid Rectangle."; Dump("", std::cerr);
  79. IsInvalid = true;
  80. }
  81. }
  82. //-----------------------------------------------------
  83. inline TVector3 ProjectionToPlane(const TVector3& point) const
  84. {
  85. Double_t d = perp * (A - point); // d = DotProduct ( N, A - point )
  86. return point + perp*d;
  87. }
  88. //-----------------------------------------------------
  89. inline LRectangle ProjectionToPlane(const LRectangle& v) const
  90. {
  91. LRectangle ret(v.volumeUID, ProjectionToPlane(v.A), ProjectionToPlane(v.B), ProjectionToPlane(v.C), ProjectionToPlane(v.D));
  92. ret. InitCenterPerp();
  93. return ret;
  94. }
  95. //-----------------------------------------------------
  96. // valid ONLY if all(point and rectangle) at the same plane
  97. inline bool IsPointInside(const TVector3& point)const
  98. {
  99. TVector3 v1 = A - B;TVector3 v2 = C - B;TVector3 v3 = point - B;
  100. double mag = v1.Mag(), proj = v3*v1/mag;
  101. if(proj < 0. || proj > mag) return false;
  102. mag = v2.Mag(), proj = v3*v2/mag;
  103. if(proj < 0. || proj > mag) return false;
  104. return true;
  105. }
  106. //-----------------------------------------------------
  107. };
  108. //------------------------------------------------------------------------------------------------------------------------
  109. class LStrip : public LRectangle
  110. {
  111. public:
  112. Int_t sectorID, detectorID, stripID;
  113. Int_t neighboring[2]; // dim same as Side_t enum
  114. TGeoCombiTrans fMatrix;
  115. LStrip();
  116. LStrip(Int_t uid, Int_t sector, Int_t detector, Int_t strip);
  117. void SetIDs(Int_t uid, Int_t sector, Int_t detector, Int_t strip){ volumeUID = uid; sectorID = sector; detectorID = detector; stripID = strip;}
  118. inline bool IsSameDetector(const LStrip& strip)const { return (sectorID == strip.sectorID && detectorID == strip.detectorID);}
  119. inline bool operator==(const LStrip& rhs)const { return (sectorID == rhs.sectorID && detectorID == rhs.detectorID && stripID == rhs.stripID);}
  120. void Dump(const char* comment = nullptr, std::ostream& out = std::cout) const;
  121. Double_t Distance(Side_t side, const LStrip& strip)const;
  122. };
  123. //------------------------------------------------------------------------------------------------------------------------
  124. class MpdTofHitProducerQA;
  125. typedef std::map<Int_t, LStrip> TmStrips; // pair<suid, LStrip>
  126. typedef Interval<LRectangle*> Tinterval;
  127. typedef IntervalTree<LRectangle*> TintervalTree;
  128. class MpdTofGeoUtils
  129. {
  130. TmStrips mStrips; //! mapping strips by suid
  131. TintervalTree mDetectorsZ, mStripsZ; //! detectors Z[cm] location interval tree
  132. TintervalTree mDetectorsPhi, mStripsPhi; //! detectors Phi[rads] location interval tree
  133. static MpdTofGeoUtils *instance;
  134. void localToMaster(const TGeoMatrix *matrix, Double_t* local, TVector3& position, std::pair<double, double>& Z, std::pair<double,double>& Phi)const;
  135. public:
  136. static MpdTofGeoUtils* Instance(){ if(instance == nullptr) instance = new MpdTofGeoUtils; return instance;}
  137. // thresh, [cm] <--- thresh. distance between neighbor strips, (see QA histo)
  138. void FindNeighborStrips(Double_t thresh, MpdTofHitProducerQA *pQA, bool forced = false);
  139. void ParseTGeoManager(MpdTofHitProducerQA* , bool forced = false, const char* flnm = nullptr);
  140. // search by strip unique ID.
  141. const LStrip* FindStrip(Int_t suid)const;
  142. bool IsPointInsideStrips(const TVector3& position, std::vector<Tinterval>& intersect, double Zerror = 0., double PhiError = 0.);
  143. bool IsPointInsideDetectors(const TVector3& position, std::vector<Tinterval>& intersect, double Zerror = 0., double PhiError = 0.);
  144. const TintervalTree* GetDetZ() { return &mDetectorsZ; };
  145. const TintervalTree* GetDetPhi() { return &mDetectorsPhi; };
  146. const TintervalTree* GetStripZ() { return &mStripsZ; };
  147. const TintervalTree* GetStripPhi() { return &mStripsPhi; };
  148. };
  149. //------------------------------------------------------------------------------------------------------------------------
  150. #endif