#include "hmdcsizescells.h"
#include "hkalgeomtools.h"
#include "TGeoTorus.h"
#include "TMath.h"
#include "TVector3.h"
#include <iostream>
using namespace std;
ClassImp(HKalGeomTools)
Double_t HKalGeomTools::distance2Points(const TVector3 &vec1, const TVector3 &vec2) {
TVector3 distV = vec1 - vec2;
return distV.Mag();
}
Double_t HKalGeomTools::distancePointToPlane(const TVector3 &point, const TVector3 &planeCenter, const TVector3 &planeNormal) {
return TMath::Abs(signedDistancePointToPlane(point, planeCenter, planeNormal));
}
Double_t HKalGeomTools::distancePointToLine(TVector3 &pca, Int_t &pcaFlag, const TVector3 &point, const TVector3 &segPoint1, const TVector3 &segPoint2) {
if(segPoint1 == segPoint2) {
::Warning("distancePointToLine()", "Both line points are equal.");
return (point - segPoint1).Mag();
}
TVector3 segDir = segPoint2 - segPoint1;
TVector3 diff = point - segPoint1;
Double_t segPar = segDir * diff / (segDir * segDir);
Double_t dist = 0.;
if(segPar > 1. || segPar <= 0.) {
pcaFlag = 1;
} else {
pcaFlag = 0;
}
pca = segPoint1 + segPar * segDir;
dist = (diff - segPar * segDir).Mag();
return dist;
}
Double_t HKalGeomTools::distancePointToSegment(const TVector3 &point, const TVector3 &segPoint1, const TVector3 &segPoint2) {
TVector3 pca;
Int_t iFlag;
return distancePointToSegment(pca, iFlag, point, segPoint1, segPoint2);
}
Double_t HKalGeomTools::distancePointToSegment(TVector3 &pca, Int_t &pcaFlag, const TVector3 &point, const TVector3 &segPoint1, const TVector3 &segPoint2) {
if(segPoint1 == segPoint2) {
::Warning("distancePointToSegment()", "Both segment end points are equal.");
return (point - segPoint1).Mag();
}
TVector3 segDir = segPoint2 - segPoint1;
TVector3 diff = point - segPoint1;
Double_t segPar = segDir * diff / (segDir * segDir);
Double_t dist = 0.;
pcaFlag = 0;
if(segPar > 0.) {
if(segPar < 1.) {
dist = (diff - segPar * segDir).Mag();
pca = segPoint1 + segPar * segDir;
} else {
dist = (diff - segDir).Mag();
pca = segPoint2;
pcaFlag = 1;
}
} else {
dist = diff.Mag();
pca = segPoint1;
pcaFlag = 1;
}
return dist;
}
Bool_t HKalGeomTools::findIntersectionLinePlane(TVector3 &pointIntersect,
const TVector3 &pos, const TVector3 &dir,
const TVector3 &planeCenter, const TVector3 &planeNormal) {
Double_t denom = planeNormal.Dot(dir);
if (denom != 0.0) {
Double_t t = ((planeCenter.x() - pos.x()) * planeNormal.x() +
(planeCenter.y() - pos.y()) * planeNormal.y() +
(planeCenter.z() - pos.z()) * planeNormal.z()) / denom;
pointIntersect = pos + (t * dir);
return kTRUE;
} else {
::Warning("findIntersection()", "No intersection point found : (plane || track)");
return kFALSE;
}
return kFALSE;
}
Bool_t HKalGeomTools::isPointOnPlane(const TVector3 &point, const TVector3 &planeCenter, const TVector3 &planeNormal){
Double_t dist = distancePointToPlane(point, planeCenter, planeNormal);
return ((TMath::Abs(dist) < 1.e-10) ? kTRUE : kFALSE);
}
Double_t HKalGeomTools::signedDistancePointToPlane(const TVector3 &point, const TVector3 &planeCenter, const TVector3 &planeNormal) {
if(planeCenter * planeNormal < 0.) {
return (point - planeCenter) * planeNormal.Unit() * -1.;
}
return (point - planeCenter) * planeNormal.Unit();
}
void HKalGeomTools::Track2ToLine(TVector3 &Pfinal, TVector3 &Pwire, Int_t &Iflag,
Double_t &Dist, Double_t &Length,
const TVector3 &X1, const TVector3 &X2,
const TVector3 &w1, const TVector3 &w2) {
TVector3 x21, x32, w21;
TVector3 xw1, xw2;
Double_t a1, b1, c1, d1, e1, t1, s1;
Double_t Delta1;
Double_t Eps = 1.E-08;
Iflag =0;
x21 = X2-X1;
w21 = w2-w1;
xw1 = X1-w1;
xw2 = X2-w1;
a1 = x21.Mag2();
b1 = x21.Dot(w21);
c1 = w21.Mag2();
d1 = xw1.Dot(x21);
e1 = xw1.Dot(w21);
Delta1 = a1*c1-b1*b1;
if(Delta1 > Eps) {
t1 = (a1*e1-b1*d1)/Delta1;
s1 = (b1*e1-c1*d1)/Delta1;
Pfinal = (X1 + x21*s1);
Pwire = (w1 + w21*t1);
Length = s1*x21.Mag();
Dist= (Pfinal-Pwire).Mag();
}
else {
Pfinal.SetXYZ(0.,0.,0.);
Pwire.SetXYZ(0.,0.,0.);
Dist=0.;
Length=0.;
Iflag = 2;
return;
}
if((((Pwire[0]<w1[0] && Pwire[0]<w2[0]) || (w2[0]<Pwire[0] && w1[0]<Pwire[0]))
&& (fabs(Pwire[0]-w1[0]) > 1e-11 && fabs(Pwire[0]- w2[0]) > 1e-11))
|| (((Pwire[1]<w1[1] && Pwire[1]<w2[1]) || (w2[1]<Pwire[1] && w1[1]<Pwire[1]))
&& (fabs(Pwire[1]-w1[1]) > 1e-11 && fabs(Pwire[1]- w2[1]) > 1e-11))
|| (((Pwire[2]<w1[2] && Pwire[2]<w2[2]) || (w2[2]<Pwire[2] && w1[2]<Pwire[2]))
&& (fabs(Pwire[2]-w1[2]) > 1e-11 && fabs(Pwire[2]- w2[2]) > 1e-11)))
{
Iflag=1;
}
}
void HKalGeomTools::TransModuleToSector(Double_t &x, Double_t &y, Double_t &z,
Int_t sector, Int_t module) {
HMdcSizesCells* fSizesCells = HMdcSizesCells::getExObject();
if(fSizesCells) {
((*fSizesCells)[sector][module]).transFrom(x,y,z);
} else {
::Error("TransModuleToSector()", "No HMdcSizesCells.");
}
}
void HKalGeomTools::TransLayerToSector(Double_t &x, Double_t &y, Double_t &z,
Int_t sector, Int_t module, Int_t layer) {
HMdcSizesCells* fSizesCells = HMdcSizesCells::getExObject();
if(fSizesCells) {
((*fSizesCells)[sector][module][layer]).transFrom(x,y,z);
} else {
::Error("TransLayerToSector()", "No HMdcSizesCells.");
}
}