/* Copyright 2008-2010, Technische Universitaet Muenchen,
Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
This file is part of GENFIT.
GENFIT is free software: you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
GENFIT is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public License
along with GENFIT. If not, see .
*/
#include "DetPlane.h"
#include
#include
#include
#include
#include
#include "TBuffer.h"
namespace genfit {
DetPlane::DetPlane(AbsFinitePlane* finite)
:finitePlane_(finite)
{
// default constructor
o_.SetXYZ(0.,0.,0.);
u_.SetXYZ(1.,0.,0.);
v_.SetXYZ(0.,1.,0.);
// sane() not needed here
}
DetPlane::DetPlane(const TVector3& o,
const TVector3& u,
const TVector3& v,
AbsFinitePlane* finite)
:o_(o), u_(u), v_(v), finitePlane_(finite)
{
sane();
}
DetPlane::DetPlane(const TVector3& o,
const TVector3& n,
AbsFinitePlane* finite)
:o_(o), finitePlane_(finite)
{
setNormal(n);
}
DetPlane::~DetPlane(){
;
}
DetPlane::DetPlane(const DetPlane& rhs) :
TObject(rhs),
o_(rhs.o_),
u_(rhs.u_),
v_(rhs.v_)
{
if(rhs.finitePlane_)
finitePlane_.reset(rhs.finitePlane_->clone());
else finitePlane_.reset();
}
DetPlane& DetPlane::operator=(DetPlane other) {
swap(other);
return *this;
}
void DetPlane::swap(DetPlane& other) {
// by swapping the members of two classes,
// the two classes are effectively swapped
std::swap(this->o_, other.o_);
std::swap(this->u_, other.u_);
std::swap(this->v_, other.v_);
this->finitePlane_.swap(other.finitePlane_);
}
void DetPlane::set(const TVector3& o,
const TVector3& u,
const TVector3& v)
{
o_ = o;
u_ = u;
v_ = v;
sane();
}
void DetPlane::setO(const TVector3& o)
{
o_ = o;
}
void DetPlane::setO(double X,double Y,double Z)
{
o_.SetXYZ(X,Y,Z);
}
void DetPlane::setU(const TVector3& u)
{
u_ = u;
sane(); // sets v_ perpendicular to u_
}
void DetPlane::setU(double X,double Y,double Z)
{
u_.SetXYZ(X,Y,Z);
sane(); // sets v_ perpendicular to u_
}
void DetPlane::setV(const TVector3& v)
{
v_ = v;
u_ = getNormal().Cross(v_);
u_ *= -1.;
sane();
}
void DetPlane::setV(double X,double Y,double Z)
{
v_.SetXYZ(X,Y,Z);
u_ = getNormal().Cross(v_);
u_ *= -1.;
sane();
}
void DetPlane::setUV(const TVector3& u,const TVector3& v)
{
u_ = u;
v_ = v;
sane();
}
void DetPlane::setON(const TVector3& o,const TVector3& n){
o_ = o;
setNormal(n);
}
TVector3 DetPlane::getNormal() const
{
return u_.Cross(v_);
}
void DetPlane::setNormal(double X,double Y,double Z){
setNormal( TVector3(X,Y,Z) );
}
void DetPlane::setNormal(const TVector3& n){
u_ = n.Orthogonal();
v_ = n.Cross(u_);
u_.SetMag(1.);
v_.SetMag(1.);
}
void DetPlane::setNormal(const double& theta, const double& phi){
setNormal( TVector3(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) );
}
TVector2 DetPlane::project(const TVector3& x)const
{
return TVector2(u_*x, v_*x);
}
TVector2 DetPlane::LabToPlane(const TVector3& x)const
{
return project(x-o_);
}
TVector3 DetPlane::toLab(const TVector2& x)const
{
TVector3 d(o_);
d += x.X()*u_;
d += x.Y()*v_;
return d;
}
TVector3 DetPlane::dist(const TVector3& x)const
{
return toLab(LabToPlane(x)) - x;
}
void DetPlane::sane(){
assert(u_!=v_);
// ensure unit vectors
u_.SetMag(1.);
v_.SetMag(1.);
// check if already orthogonal
if (u_.Dot(v_) < 1.E-5) return;
// ensure orthogonal system
v_ = getNormal().Cross(u_);
}
void DetPlane::Print(const Option_t* option) const
{
std::cout<<"DetPlane: "
<<"O("<Print(option);
}
}
/*
I could write pages of comments about correct equality checking for
floating point numbers, but: When two planes are as close as 10E-5 cm
in all nine numbers that define the plane, this will be enough for all
practical purposes
*/
bool operator== (const DetPlane& lhs, const DetPlane& rhs){
if (&lhs == &rhs)
return true;
static const double detplaneEpsilon = 1.E-5;
if(
fabs( (lhs.o_.X()-rhs.o_.X()) ) > detplaneEpsilon ||
fabs( (lhs.o_.Y()-rhs.o_.Y()) ) > detplaneEpsilon ||
fabs( (lhs.o_.Z()-rhs.o_.Z()) ) > detplaneEpsilon
) return false;
else if(
fabs( (lhs.u_.X()-rhs.u_.X()) ) > detplaneEpsilon ||
fabs( (lhs.u_.Y()-rhs.u_.Y()) ) > detplaneEpsilon ||
fabs( (lhs.u_.Z()-rhs.u_.Z()) ) > detplaneEpsilon
) return false;
else if(
fabs( (lhs.v_.X()-rhs.v_.X()) ) > detplaneEpsilon ||
fabs( (lhs.v_.Y()-rhs.v_.Y()) ) > detplaneEpsilon ||
fabs( (lhs.v_.Z()-rhs.v_.Z()) ) > detplaneEpsilon
) return false;
return true;
}
bool operator!= (const DetPlane& lhs, const DetPlane& rhs){
return !(lhs==rhs);
}
double DetPlane::distance(const TVector3& point) const {
// |(point - o_)*(u_ x v_)|
return fabs( (point.X()-o_.X()) * (u_.Y()*v_.Z() - u_.Z()*v_.Y()) +
(point.Y()-o_.Y()) * (u_.Z()*v_.X() - u_.X()*v_.Z()) +
(point.Z()-o_.Z()) * (u_.X()*v_.Y() - u_.Y()*v_.X()));
}
double DetPlane::distance(double x, double y, double z) const {
// |(point - o_)*(u_ x v_)|
return fabs( (x-o_.X()) * (u_.Y()*v_.Z() - u_.Z()*v_.Y()) +
(y-o_.Y()) * (u_.Z()*v_.X() - u_.X()*v_.Z()) +
(z-o_.Z()) * (u_.X()*v_.Y() - u_.Y()*v_.X()));
}
TVector2 DetPlane::straightLineToPlane (const TVector3& point, const TVector3& dir) const {
TVector3 dirNorm(dir.Unit());
TVector3 normal = getNormal();
double dirTimesN = dirNorm*normal;
if(fabs(dirTimesN)<1.E-6){//straight line is parallel to plane, so return infinity
return TVector2(1.E100,1.E100);
}
double t = 1./dirTimesN * ((o_-point)*normal);
return project(point - o_ + t * dirNorm);
}
//! gives u,v coordinates of the intersection point of a straight line with plane
void DetPlane::straightLineToPlane(const double& posX, const double& posY, const double& posZ,
const double& dirX, const double& dirY, const double& dirZ,
double& u, double& v) const {
TVector3 W = getNormal();
double dirTimesN = dirX*W.X() + dirY*W.Y() + dirZ*W.Z();
if(fabs(dirTimesN)<1.E-6){//straight line is parallel to plane, so return infinity
u = 1.E100;
v = 1.E100;
return;
}
double t = 1./dirTimesN * ((o_.X()-posX)*W.X() +
(o_.Y()-posY)*W.Y() +
(o_.Z()-posZ)*W.Z());
double posOnPlaneX = posX-o_.X() + t*dirX;
double posOnPlaneY = posY-o_.Y() + t*dirY;
double posOnPlaneZ = posZ-o_.Z() + t*dirZ;
u = u_.X()*posOnPlaneX + u_.Y()*posOnPlaneY + u_.Z()*posOnPlaneZ;
v = v_.X()*posOnPlaneX + v_.Y()*posOnPlaneY + v_.Z()*posOnPlaneZ;
}
void DetPlane::rotate(double angle) {
TVector3 normal = getNormal();
u_.Rotate(angle, normal);
v_.Rotate(angle, normal);
sane();
}
void DetPlane::reset() {
o_.SetXYZ(0.,0.,0.);
u_.SetXYZ(1.,0.,0.);
v_.SetXYZ(0.,1.,0.);
finitePlane_.reset();
}
void DetPlane::Streamer(TBuffer &R__b)
{
// Stream an object of class genfit::DetPlane.
// This is modified from the streamer generated by ROOT in order to
// account for the scoped_ptr.
//This works around a msvc bug and should be harmless on other platforms
typedef ::genfit::DetPlane thisClass;
UInt_t R__s, R__c;
if (R__b.IsReading()) {
Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
//TObject::Streamer(R__b);
o_.Streamer(R__b);
u_.Streamer(R__b);
v_.Streamer(R__b);
finitePlane_.reset();
char flag;
R__b >> flag;
if (flag) {
// Read AbsFinitePlane with correct overload ... see comment
// in writing code.
TClass* cl = TClass::Load(R__b);
AbsFinitePlane *p = (AbsFinitePlane*)(cl->New());
cl->Streamer(p, R__b);
finitePlane_.reset(p);
}
R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
} else {
R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
//TObject::Streamer(R__b);
o_.Streamer(R__b);
u_.Streamer(R__b);
v_.Streamer(R__b);
if (finitePlane_) {
R__b << (char)1;
// finitPlane_ is a scoped_ptr, but a shared plane may be
// written several times in different places (e.g. it may also
// be stored in the measurement class). Since we have no way
// of knowing in which other places the same plane is written
// (it may even be in another track!), we cannot synchronize
// these writes and have to duplicate the SharedPlanePtr's
// instead. ROOT caches which pointers were written and read
// if one uses 'R__b << p' or equivalent. This can lead to
// trouble have no way of synchronizing two reads to
// shared_ptr's pointing to the same memory location. But
// even if we break the link between the two shared_ptr's,
// ROOT will still try to outsmart us, and it will notice that
// the finitePlane_ was the same during writing. This will
// lead to the same address being attached to two different
// scoped_ptrs in reading. Highly undesirable. Since we
// cannot know whether other parts of code implicitly or
// explicitly rely on TBuffer's caching, we cannot just
// disable this caching via R__b.ResetMap() (it's not
// documented in any elucidating manner anyway).
//
// Therefore, we have to write and read the AbsFinitePlane*
// manually, bypassing ROOT's caching. In order to do this,
// we need the information on the actual type, because
// otherwise we couldn't read it back reliably. Finally, the
// _working_ means of reading and writing class information
// are TClass::Store and TClass::Load, again barely
// documented, but if one tries any of the other ways that
// suggest themselves, weird breakage will occur.
finitePlane_->IsA()->Store(R__b);
finitePlane_->Streamer(R__b);
} else {
R__b << (char)0;
}
R__b.SetByteCount(R__c, kTRUE);
}
}
} /* End of namespace genfit */