sPhenix code displayed by LXR master/​GenFit/​core/​src/​DetPlane.cc	Source navigation Diff markup Identifier search General search
	Version: master Revert   Change

File indexing completed on 2025-08-05 08:18:21

 /* 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 <http://www.gnu.org/licenses/>.
 */
 
 #include "DetPlane.h"
 #include "IO.h"
 
 #include <cassert>
 #include <cmath>
 #include <TMath.h>
 #include <TClass.h>
 #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
 {
   printOut<<"DetPlane: "
      <<"O("<<o_.X()<<", "<<o_.Y()<<", "<<o_.Z()<<") "
      <<"u("<<u_.X()<<", "<<u_.Y()<<", "<<u_.Z()<<") "
      <<"v("<<v_.X()<<", "<<v_.Y()<<", "<<v_.Z()<<") "
      <<"n("<<getNormal().X()<<", "<<getNormal().Y()<<", "<<getNormal().Z()<<") "
        <<std::endl;
   if(finitePlane_ != nullptr) {
     finitePlane_->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 */

This page was automatically generated by the 2.3.7 LXR engine. The LXR team