MpdTofGeoUtils.cxx 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512
  1. //------------------------------------------------------------------------------------------------------------------------
  2. #include <assert.h>
  3. #include <iostream>
  4. #include <fstream>
  5. #include <TGeoManager.h>
  6. #include <TRandom2.h>
  7. #include <TRotation.h>
  8. #include <TClonesArray.h>
  9. #include <TH1D.h>
  10. #include <TH2D.h>
  11. #include <TMath.h>
  12. #include <TGeoMatrix.h>
  13. #include <TGeoBBox.h>
  14. #include <TGeoNode.h>
  15. #include "FairRootManager.h"
  16. #include "FairLogger.h"
  17. #include "MpdTof.h"
  18. #include "MpdTofHitProducerQA.h"
  19. #include "MpdTofPoint.h"
  20. #include "MpdTofGeoUtils.h"
  21. using namespace std;
  22. MpdTofGeoUtils* MpdTofGeoUtils::instance = nullptr;
  23. //------------------------------------------------------------------------------------------------------------------------
  24. void MpdTofGeoUtils::FindNeighborStrips(Double_t thresh, MpdTofHitProducerQA *pQA, bool forced)
  25. {
  26. static bool found = false;
  27. if( !forced && found) return; // neighbor strips already found and filled
  28. size_t NR = 0, NL= 0;
  29. for(auto& it1 : mStrips) // cycle1 by strips
  30. {
  31. LStrip *strip1 = &(it1.second);
  32. for(const auto& it2 : mStrips) // cycle2 by strips
  33. {
  34. const LStrip *strip2 = &(it2.second);
  35. if(strip1 == strip2) continue; // skip calculate yourself
  36. // CATION: Ckeck only left and right sides(one row strips NOW)
  37. double distance = strip1->Distance(LStrip::kRight, *strip2);
  38. if(pQA) pQA->FillDistanceToStrip(distance);
  39. if(distance < thresh && MpdTofPoint::IsSameGap(strip1->volumeUID, strip2->volumeUID))
  40. {
  41. strip1->neighboring[LStrip::kRight] = strip2->volumeUID; NR++;
  42. if(pQA) pQA->FillNeighborPairs(strip2->stripID, strip1->stripID);
  43. }
  44. distance = strip1->Distance(LStrip::kLeft, *strip2);
  45. if(pQA) pQA->FillDistanceToStrip(distance* -1.);
  46. if(distance < thresh && MpdTofPoint::IsSameGap(strip1->volumeUID, strip2->volumeUID))
  47. {
  48. strip1->neighboring[LStrip::kLeft] = strip2->volumeUID; NL++;
  49. if(pQA) pQA->FillNeighborPairs(strip2->stripID, strip1->stripID);
  50. }
  51. } // cycle2 by strips
  52. } // cycle1 by strips
  53. found = true;
  54. LOG(DEBUG1)<<"[MpdTofGeoUtils::FindNeighborStrips] Neighbor strips: left = "<<NL<<", right = "<<NR<<".";
  55. }
  56. //------------------------------------------------------------------------------------------------------------------------
  57. void MpdTofGeoUtils::ParseTGeoManager(MpdTofHitProducerQA *pQA, bool forced, const char* flnm)
  58. {
  59. constexpr double degree180 = TMath::Pi();
  60. assert(gGeoManager);
  61. if( !forced && !mStrips.empty()) return; // already parsed and filled
  62. mStrips.clear();
  63. ofstream ofs;
  64. if(flnm) ofs.open(flnm);
  65. const TString pathTOF = "/cave_1/tof1_0";
  66. gGeoManager->cd(pathTOF);
  67. Double_t *local = new Double_t[3](), master[3] = {0}, dX, dY, dZ;
  68. int nSectors = 0, nDetectors = 0, nStrips = 0;
  69. TObjArray *array = gGeoManager->GetCurrentVolume()->GetNodes();
  70. auto it1 = array->MakeIterator();
  71. vector<Tinterval> vDetectorsZ, vDetectorsPhi;
  72. vector<Tinterval> vStripsZ, vStripsPhi;
  73. while( TGeoNode *sectorNode = (TGeoNode*) it1->Next() ) // SECTORS vSector_#, DetectorBoxV_# +, StripBoxV_#, StripActiveGasV_# +
  74. {
  75. if(!TString(sectorNode->GetName()).Contains("tof1Sector")) continue;
  76. TString PATH1 = pathTOF + "/" + sectorNode->GetName();
  77. Int_t sectorID = sectorNode->GetNumber(); // sector [1,...,14]
  78. nSectors++;
  79. if(ofs.is_open()) ofs<<"\n SECTOR: "<<sectorNode->GetName()<<", copy# "<<sectorID<<" path= "<<PATH1.Data();
  80. auto it3 = sectorNode->GetNodes()->MakeIterator();
  81. while( TGeoNode *detectorNode = (TGeoNode*) it3->Next() ) // DETECTORS
  82. {
  83. TString detName = detectorNode->GetName();
  84. if(!detName.Contains("DetectorBox")) continue;
  85. TString PATH3 = PATH1 + "/" + detName;
  86. Int_t detectorID = detectorNode->GetNumber(); // detector [1,...,20]
  87. nDetectors++;
  88. if(ofs.is_open()) ofs<<"\n\t\t DETECTOR: "<<detName.Data()<<", copy# "<<detectorID<<" path= "<<PATH3.Data();
  89. gGeoManager->cd(PATH3); // cd to detector box node
  90. TGeoMatrix *matrix = gGeoManager->GetCurrentMatrix();
  91. TGeoBBox *box = (TGeoBBox*) gGeoManager->GetCurrentNode()->GetVolume()->GetShape();
  92. dX = box->GetDX(); dY = box->GetDY(); dZ = box->GetDZ();
  93. TVector3 A,B,C,D; // detector edges
  94. auto Z_range = make_pair(1.e+10, -1.e+10), Phi_range(Z_range);
  95. local[0] = -dX; local[1] = -dY +1.3-0.26; local[2] = +dZ; // shift by Y to 1th gas sensitive layer of strips (for rough track propagate to plate)
  96. localToMaster(matrix, local, A, Z_range, Phi_range);
  97. local[0] = +dX; local[1] = -dY +1.3-0.26; local[2] = +dZ;
  98. localToMaster(matrix, local, B, Z_range, Phi_range);
  99. local[0] = +dX; local[1] = -dY +1.3-0.26; local[2] = -dZ;
  100. localToMaster(matrix, local, C, Z_range, Phi_range);
  101. local[0] = -dX; local[1] = -dY +1.3-0.26; local[2] = -dZ;
  102. localToMaster(matrix, local, D, Z_range, Phi_range);
  103. if(pQA) pQA->FillDetectorsMap(A, B, C, D);
  104. // Add Detector location interval to intervalTree
  105. Int_t duid = MpdTofPoint::ClearGap(MpdTofPoint::GetSuid72(sectorID, detectorID, 0)); // suid for 0th strip of the selected detector, gap reset to 0
  106. LRectangle * det = new LRectangle(duid, A, B, C, D);
  107. det->InitCenterPerp();
  108. vDetectorsZ.push_back(Tinterval(Z_range.first, Z_range.second, det));
  109. if(Phi_range.second - Phi_range.first > degree180) // segment overlap the boundary [-180.,+180[
  110. {
  111. vDetectorsPhi.push_back(Tinterval(-degree180, Phi_range.first, det));
  112. vDetectorsPhi.push_back(Tinterval(Phi_range.second, degree180, det));
  113. cout<<"\n OVERLAPED phimin="<<Phi_range.first<<", phimax="<<Phi_range.second;
  114. }
  115. else vDetectorsPhi.push_back(Tinterval(Phi_range.first, Phi_range.second, det));
  116. auto it4 = detectorNode->GetNodes()->MakeIterator();
  117. while( TGeoNode *stripBoxNode = (TGeoNode*) it4->Next() ) // STRIPS
  118. {
  119. TString stripBoxName = stripBoxNode->GetName();
  120. if(!stripBoxName.Contains("StripActiveGas")) continue;
  121. TString PATH4 = PATH3 + "/" + stripBoxName;
  122. Int_t stripID = stripBoxNode->GetNumber(); // strip[1,...,72]
  123. nStrips++;
  124. if(ofs.is_open()) ofs<<"\n\t\t\t STRIP: "<<stripBoxName.Data()<<", copy# "<<stripID<<" path= "<<PATH4.Data();
  125. gGeoManager->cd(PATH4); // cd to strip sensitive gas node
  126. TGeoMatrix *matrix4 = gGeoManager->GetCurrentMatrix();
  127. TGeoBBox *box4 = (TGeoBBox*) gGeoManager->GetCurrentNode()->GetVolume()->GetShape();
  128. dX = box4->GetDX(); dY = box4->GetDY(); dZ = box4->GetDZ();
  129. Int_t suid = MpdTofPoint::GetSuid72(sectorID, detectorID, stripID);
  130. if(ofs.is_open())
  131. {
  132. MpdTofPoint::PrintSuid(suid, "\n\t\t\t\t Add strip:", ofs);
  133. ofs<<" center= "<<master[0]<<", "<<master[1]<<", "<<master[2]<<", dXYZ= "<<dX<<", "<<dY<<", "<<dZ;
  134. }
  135. LStrip strip(suid, sectorID, detectorID, stripID);
  136. strip.fMatrix = (*matrix4);
  137. auto sZ_range = make_pair(1.e+10, -1.e+10), sPhi_range(sZ_range); // init by big values
  138. //
  139. // ^ +Z
  140. // | A right side B
  141. // | |----------------------------------|
  142. // | |-------------strip----------------|
  143. // | D left side C
  144. // |-------------------------------------------------> +X
  145. // | -Z
  146. // edges on the front plate of the strips. perp along Y.
  147. local[0] = -dX; local[1] = -dY; local[2] = +dZ;
  148. localToMaster(matrix4, local, strip.A, sZ_range, sPhi_range);
  149. local[0] = +dX; local[1] = -dY; local[2] = +dZ;
  150. localToMaster(matrix4, local, strip.B, sZ_range, sPhi_range);
  151. local[0] = +dX; local[1] = -dY; local[2] = -dZ;
  152. localToMaster(matrix4, local, strip.C, sZ_range, sPhi_range);
  153. local[0] = -dX; local[1] = -dY; local[2] = -dZ;
  154. localToMaster(matrix4, local, strip.D, sZ_range, sPhi_range);
  155. strip.InitCenterPerp(); // calc. strip center and perp. AFTER setup A,B,C,D
  156. if(pQA) pQA->FillStripsMap(strip.A, strip.B, strip.C, strip.D);
  157. auto pStrip = new LStrip(strip);
  158. vStripsZ.push_back(Tinterval(sZ_range.first, sZ_range.second, pStrip));
  159. if(sPhi_range.second - sPhi_range.first > degree180) // segment overlap the boundary [-180.,+180[
  160. {
  161. vStripsPhi.push_back(Tinterval(-degree180, sPhi_range.first, pStrip));
  162. vStripsPhi.push_back(Tinterval(sPhi_range.second, degree180, pStrip));
  163. if(ofs.is_open()) ofs<<"\n OVERLAPED strip phimin="<<sPhi_range.first<<", phimax="<<sPhi_range.second;
  164. }
  165. else vStripsPhi.push_back(Tinterval(sPhi_range.first, sPhi_range.second, pStrip));
  166. ///strip.Dump("\n strip:");
  167. bool IsUniqueUID = mStrips.insert(make_pair(suid, strip)).second;
  168. assert(IsUniqueUID);
  169. } // STRIPS
  170. } // DETECTORS
  171. } // SECTORS
  172. mStripsZ = TintervalTree(vStripsZ);
  173. mStripsPhi = TintervalTree(vStripsPhi);
  174. mDetectorsZ = TintervalTree(vDetectorsZ);
  175. mDetectorsPhi = TintervalTree(vDetectorsPhi);
  176. if(ofs.is_open()) ofs.close();
  177. delete[] local;
  178. LOG(DEBUG1)<<"[MpdTof::ParseTGeoManager] Sectors="<<nSectors<<", detectors="<<nDetectors<<", strips="<<nStrips<<".";
  179. }
  180. //------------------------------------------------------------------------------------------------------------------------
  181. void MpdTofGeoUtils::localToMaster(const TGeoMatrix *matrix, Double_t* local, TVector3& position, pair<double, double>& Z, pair<double,double>& Phi)const
  182. {
  183. Double_t master[3];
  184. matrix->LocalToMaster(local, master);
  185. position.SetXYZ(master[0], master[1], master[2]);
  186. Z.first = std::min(position.Z(), Z.first);
  187. Z.second = std::max(position.Z(), Z.second);
  188. Phi.first = std::min(position.Phi(), Phi.first);
  189. Phi.second = std::max(position.Phi(), Phi.second);
  190. }
  191. //------------------------------------------------------------------------------------------------------------------------
  192. const LStrip* MpdTofGeoUtils::FindStrip(Int_t suid) const
  193. {
  194. auto cit = mStrips.find(suid);
  195. assert(!(cit == mStrips.end() && MpdTofPoint::PrintSuid(suid, "\n [MpdTofGeoUtils::FindStrip] : unknown suid. ", cerr)));
  196. return &(cit->second);
  197. }
  198. //------------------------------------------------------------------------------------------------------------------------
  199. bool MpdTofGeoUtils::IsPointInsideStrips(const TVector3& position, vector<Tinterval>& intersect, double Zerror, double PhiError)
  200. {
  201. const double Z = position.Z();
  202. const double Phi = position.Phi();
  203. vector<Tinterval> insideZ, insidePhi;
  204. mStripsZ.findOverlapping(Z - Zerror, Z + Zerror, insideZ);
  205. mStripsPhi.findOverlapping(Phi - PhiError, Phi + PhiError, insidePhi);
  206. if(!insideZ.empty() && !insidePhi.empty()) // have overlaped segments both Z and Phi
  207. {
  208. intersect.clear();
  209. // calc. intersection
  210. sort(insideZ.begin(), insideZ.end(), IntervalValueSorter<LRectangle*>());
  211. sort(insidePhi.begin(), insidePhi.end(), IntervalValueSorter<LRectangle*>());
  212. set_intersection(insideZ.begin(), insideZ.end(), insidePhi.begin(), insidePhi.end(), back_inserter(intersect), IntervalValueSorter<LRectangle*>());
  213. return true;
  214. }
  215. return false;
  216. }
  217. //------------------------------------------------------------------------------------------------------------------------
  218. bool MpdTofGeoUtils::IsPointInsideDetectors(const TVector3& position, vector<Tinterval>& intersect, double Zerror, double PhiError)
  219. {
  220. const double Z = position.Z();
  221. const double Phi = position.Phi();
  222. vector<Tinterval> insideZ, insidePhi;
  223. mDetectorsZ.findOverlapping(Z - Zerror, Z + Zerror, insideZ);
  224. mDetectorsPhi.findOverlapping(Phi - PhiError, Phi + PhiError, insidePhi);
  225. if(!insideZ.empty() && !insidePhi.empty()) // have overlaped segments both Z and Phi
  226. {
  227. intersect.clear();
  228. // calc. intersection
  229. sort(insideZ.begin(), insideZ.end(), IntervalValueSorter<LRectangle*>());
  230. sort(insidePhi.begin(), insidePhi.end(), IntervalValueSorter<LRectangle*>());
  231. set_intersection(insideZ.begin(), insideZ.end(), insidePhi.begin(), insidePhi.end(), back_inserter(intersect), IntervalValueSorter<LRectangle*>());
  232. return true;
  233. }
  234. return false;
  235. }
  236. //------------------------------------------------------------------------------------------------------------------------
  237. //------------------------------------------------------------------------------------------------------------------------
  238. //------------------------------------------------------------------------------------------------------------------------
  239. LRectangle::LRectangle(Int_t uid, const TVector3& a, const TVector3& b, const TVector3& c, const TVector3& d, bool check)
  240. : IsInvalid(false), volumeUID(uid), A(a), B(b), C(c), D(d)
  241. {
  242. if(check) CheckInValid();
  243. }
  244. //------------------------------------------------------------------------------------------------------------------------
  245. void LRectangle::GetRPhiRanges(Double_t& Rmin, Double_t& Rmax, Double_t& PhiMin, Double_t& PhiMax)
  246. {
  247. Rmin = PhiMin = 1.e+10;
  248. Rmax = PhiMax = -1.e+10;
  249. getRPhi(A, Rmin, Rmax, PhiMin, PhiMax); // edge
  250. getRPhi(B, Rmin, Rmax, PhiMin, PhiMax);
  251. getRPhi(C, Rmin, Rmax, PhiMin, PhiMax);
  252. getRPhi(D, Rmin, Rmax, PhiMin, PhiMax);
  253. getRPhi((A+B)*0.5, Rmin, Rmax, PhiMin, PhiMax); // middle of side
  254. getRPhi((B+C)*0.5, Rmin, Rmax, PhiMin, PhiMax);
  255. getRPhi((C+D)*0.5, Rmin, Rmax, PhiMin, PhiMax);
  256. getRPhi((D+A)*0.5, Rmin, Rmax, PhiMin, PhiMax);
  257. }
  258. //------------------------------------------------------------------------------------------------------------------------
  259. Double_t LRectangle::DistanceFromPointToLine(const TVector3* pos, const TVector3& P1,const TVector3& P2)const
  260. {
  261. assert(P1 != P2);
  262. return ( (*pos - P1).Cross(*pos - P2) ).Mag() / (P2 - P1).Mag();
  263. }
  264. //------------------------------------------------------------------------------------------------------------------------
  265. void LRectangle::Rotate(const TRotation& rot)
  266. {
  267. A *= rot;
  268. B *= rot;
  269. C *= rot;
  270. D *= rot;
  271. center *= rot;
  272. perp *= rot;
  273. }
  274. //------------------------------------------------------------------------------------------------------------------------
  275. Double_t LRectangle::DistanceFromPointToLineSegment(const TVector3* pos, const TVector3& P1,const TVector3& P2)const
  276. {
  277. assert(P1 != P2);
  278. TVector3 v = P2 - P1;
  279. TVector3 w = (*pos) - P1;
  280. double c1 = w.Dot(v);
  281. if( c1 <= 0 ) return w.Mag();
  282. double c2 = v.Dot(v);
  283. if( c2 <= c1 ) return ((*pos) - P2).Mag();
  284. TVector3 Pb = P1 + (c1/c2) * v;
  285. return ((*pos) - Pb).Mag();
  286. }
  287. //------------------------------------------------------------------------------------------------------------------------
  288. Double_t LRectangle::MinDistanceToEdge(const TVector3* pos, Side_t& side) const
  289. {
  290. double right = DistanceFromPointToLineSegment(pos, A, B);
  291. double left = DistanceFromPointToLineSegment(pos, C, D);
  292. // sorting & return minimal value
  293. if( right <= left )
  294. {
  295. side = LStrip::kRight;
  296. return right;
  297. }
  298. side = LStrip::kLeft;
  299. return left;
  300. }
  301. //------------------------------------------------------------------------------------------------------------------------
  302. void LRectangle::Dump(const char* comment, ostream& out) const
  303. {
  304. if(comment) out<<comment;
  305. out<<" uid="<<volumeUID<<" IsInvalid="<<IsInvalid;
  306. MpdTof::Print(A, " A:", out); MpdTof::Print(B, " B:", out); MpdTof::Print(C, " C:", out); MpdTof::Print(D, " D:", out); MpdTof::Print(perp, " perp:", out);
  307. }
  308. //------------------------------------------------------------------------------------------------------------------------
  309. void LRectangle::Test(const char* comment, ostream &out)
  310. {
  311. if(comment) out<<comment;
  312. out<<"\n [LRectangle::Test]--------------------------->>>";
  313. //const double epsilon = 1.E-3;
  314. for(double xshift = -10; xshift < 10; xshift++)
  315. for(double yshift = -10; yshift < 10; yshift++)
  316. {
  317. LRectangle A(1, TVector3(0,0,0), TVector3(0,10,0), TVector3(10,10,0), TVector3(10,0,0)); A.InitCenterPerp(); // 10x10
  318. LRectangle B(2, TVector3(0,0,1), TVector3(0,5,1), TVector3(5,5,1), TVector3(5,0,1)); B.InitCenterPerp(); // 5x5
  319. B.Shift({xshift, yshift, 0.});
  320. double square = B.GetIntersectionSquare(A);// smaller intersect biger
  321. out<<"\n shift: ("<<xshift<<","<<yshift<<") square="<<square;
  322. // Randomize rectangle orientation
  323. TRotation rot; TVector3 axis, shift;
  324. for(auto i=0;i<100;i++)
  325. {
  326. axis.SetXYZ(gRandom->Uniform(), gRandom->Uniform(), gRandom->Uniform());
  327. shift.SetXYZ(gRandom->Uniform(-10,10), gRandom->Uniform(-10,10), gRandom->Uniform(-10,10));
  328. rot.Rotate(gRandom->Uniform(-10.,10.), axis); // random rotate around axis
  329. A.Rotate(rot);
  330. B.Rotate(rot);
  331. A.Shift(shift);
  332. B.Shift(shift);
  333. double newsquare = B.GetIntersectionSquare(A);
  334. if(fabs(newsquare - square) > 1.) out<<"\n\tperp.: theta="<<A.perp.Theta()*TMath::RadToDeg()<<" phi="<<A.perp.Phi()*TMath::RadToDeg()<<", int. square="<<newsquare;
  335. }
  336. }
  337. }
  338. //------------------------------------------------------------------------------------------------------------------------
  339. Double_t LRectangle::GetIntersectionSquare(const LRectangle& v1)const
  340. {
  341. LRectangle v(v1);
  342. if(!IsSamePlane(v1).first)
  343. {
  344. v = ProjectionToPlane(v1);
  345. }
  346. if(!IsSamePlane(v).first) return 0.;
  347. TVector3 p, AB = B - A, AD = D - A;
  348. const size_t Nprobes = 10000;
  349. size_t Ncounter = 0;
  350. for(size_t i=0;i<Nprobes;i++)
  351. {
  352. p = A + AB * gRandom->Uniform() + AD * gRandom->Uniform();
  353. if(v.IsPointInside(p)) Ncounter++;
  354. }
  355. if(Ncounter == 0) return 0.;
  356. return Square()*Ncounter/((double)Nprobes);
  357. }
  358. //------------------------------------------------------------------------------------------------------------------------
  359. //------------------------------------------------------------------------------------------------------------------------
  360. //------------------------------------------------------------------------------------------------------------------------
  361. LStrip::LStrip()
  362. : LRectangle(), sectorID(kInvalid), detectorID(kInvalid), stripID(kInvalid)
  363. {
  364. neighboring[kRight] = kInvalid;
  365. neighboring[kLeft] = kInvalid;
  366. }
  367. //------------------------------------------------------------------------------------------------------------------------
  368. LStrip::LStrip(Int_t uid, Int_t sector, Int_t detector, Int_t strip)
  369. : LRectangle(), sectorID(sector), detectorID(detector), stripID(strip)
  370. {
  371. volumeUID = uid;
  372. neighboring[kRight] = kInvalid;
  373. neighboring[kLeft] = kInvalid;
  374. }
  375. //------------------------------------------------------------------------------------------------------------------------
  376. void LStrip::Dump(const char* comment, ostream& out) const
  377. {
  378. if(comment) out<<comment;
  379. out<<" ids: "<<sectorID<<", "<<detectorID<<", "<<stripID;
  380. LRectangle::Dump(nullptr, out);
  381. }
  382. //------------------------------------------------------------------------------------------------------------------------
  383. Double_t LStrip::Distance(Side_t side, const LStrip& strip) const
  384. {
  385. Double_t value, min1 = 1.e+10, min2 = 1.e+10; // big value
  386. if((*this) == strip) return min1; // same strip, return big value
  387. if(!IsSameDetector(strip)) return min1; // different detector, return big value
  388. const TVector3 *p1{nullptr}, *p2{nullptr}; //TODO test
  389. switch(side)
  390. {
  391. case kRight: p1 = &A; p2 = &B; break;
  392. case kLeft: p1 = &C; p2 = &D; break;
  393. case kInvalid: return min1;
  394. }
  395. value = (*p1 - strip.A).Mag(); min1 = std::min(value, min1);
  396. value = (*p1 - strip.B).Mag(); min1 = std::min(value, min1);
  397. value = (*p1 - strip.C).Mag(); min1 = std::min(value, min1);
  398. value = (*p1 - strip.D).Mag(); min1 = std::min(value, min1);
  399. value = (*p2 - strip.A).Mag(); min2 = std::min(value, min2);
  400. value = (*p2 - strip.B).Mag(); min2 = std::min(value, min2);
  401. value = (*p2 - strip.C).Mag(); min2 = std::min(value, min2);
  402. value = (*p2 - strip.D).Mag(); min2 = std::min(value, min2);
  403. return std::min(min1, min2);
  404. }
  405. //------------------------------------------------------------------------------------------------------------------------