Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DetPlane.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DetPlane.cc
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "DetPlane.h"
21 #include "IO.h"
22 
23 #include <cassert>
24 #include <cmath>
25 #include <TMath.h>
26 #include <TClass.h>
27 #include <TBuffer.h>
28 
29 namespace genfit {
30 
31 
33  :finitePlane_(finite)
34 {
35  // default constructor
36  o_.SetXYZ(0.,0.,0.);
37  u_.SetXYZ(1.,0.,0.);
38  v_.SetXYZ(0.,1.,0.);
39  // sane() not needed here
40 }
41 
42 
43 DetPlane::DetPlane(const TVector3& o,
44  const TVector3& u,
45  const TVector3& v,
46  AbsFinitePlane* finite)
47  :o_(o), u_(u), v_(v), finitePlane_(finite)
48 {
49  sane();
50 }
51 
52 DetPlane::DetPlane(const TVector3& o,
53  const TVector3& n,
54  AbsFinitePlane* finite)
55  :o_(o), finitePlane_(finite)
56 {
57  setNormal(n);
58 }
59 
60 
62  ;
63 }
64 
65 
67  TObject(rhs),
68  o_(rhs.o_),
69  u_(rhs.u_),
70  v_(rhs.v_)
71 {
72  if(rhs.finitePlane_)
73  finitePlane_.reset(rhs.finitePlane_->clone());
74  else finitePlane_.reset();
75 }
76 
77 
79  swap(other);
80  return *this;
81 }
82 
83 
84 void DetPlane::swap(DetPlane& other) {
85  // by swapping the members of two classes,
86  // the two classes are effectively swapped
87  std::swap(this->o_, other.o_);
88  std::swap(this->u_, other.u_);
89  std::swap(this->v_, other.v_);
90  this->finitePlane_.swap(other.finitePlane_);
91 }
92 
93 
94 void DetPlane::set(const TVector3& o,
95  const TVector3& u,
96  const TVector3& v)
97 {
98  o_ = o;
99  u_ = u;
100  v_ = v;
101  sane();
102 }
103 
104 
105 void DetPlane::setO(const TVector3& o)
106 {
107  o_ = o;
108 }
109 
110 void DetPlane::setO(double X,double Y,double Z)
111 {
112  o_.SetXYZ(X,Y,Z);
113 }
114 
115 void DetPlane::setU(const TVector3& u)
116 {
117  u_ = u;
118  sane(); // sets v_ perpendicular to u_
119 }
120 
121 void DetPlane::setU(double X,double Y,double Z)
122 {
123  u_.SetXYZ(X,Y,Z);
124  sane(); // sets v_ perpendicular to u_
125 }
126 
127 void DetPlane::setV(const TVector3& v)
128 {
129  v_ = v;
130  u_ = getNormal().Cross(v_);
131  u_ *= -1.;
132  sane();
133 }
134 
135 void DetPlane::setV(double X,double Y,double Z)
136 {
137  v_.SetXYZ(X,Y,Z);
138  u_ = getNormal().Cross(v_);
139  u_ *= -1.;
140  sane();
141 }
142 
143 void DetPlane::setUV(const TVector3& u,const TVector3& v)
144 {
145  u_ = u;
146  v_ = v;
147  sane();
148 }
149 
150 void DetPlane::setON(const TVector3& o,const TVector3& n){
151  o_ = o;
152  setNormal(n);
153 }
154 
155 
156 TVector3 DetPlane::getNormal() const
157 {
158  return u_.Cross(v_);
159 }
160 
161 void DetPlane::setNormal(double X,double Y,double Z){
162  setNormal( TVector3(X,Y,Z) );
163 }
164 
165 void DetPlane::setNormal(const TVector3& n){
166  u_ = n.Orthogonal();
167  v_ = n.Cross(u_);
168  u_.SetMag(1.);
169  v_.SetMag(1.);
170 }
171 
172 void DetPlane::setNormal(const double& theta, const double& phi){
173  setNormal( TVector3(TMath::Sin(theta)*TMath::Cos(phi),TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta)) );
174 }
175 
176 
177 TVector2 DetPlane::project(const TVector3& x)const
178 {
179  return TVector2(u_*x, v_*x);
180 }
181 
182 
183 TVector2 DetPlane::LabToPlane(const TVector3& x)const
184 {
185  return project(x-o_);
186 }
187 
188 
189 TVector3 DetPlane::toLab(const TVector2& x)const
190 {
191  TVector3 d(o_);
192  d += x.X()*u_;
193  d += x.Y()*v_;
194  return d;
195 }
196 
197 
198 TVector3 DetPlane::dist(const TVector3& x)const
199 {
200  return toLab(LabToPlane(x)) - x;
201 }
202 
203 
205  assert(u_!=v_);
206 
207  // ensure unit vectors
208  u_.SetMag(1.);
209  v_.SetMag(1.);
210 
211  // check if already orthogonal
212  if (u_.Dot(v_) < 1.E-5) return;
213 
214  // ensure orthogonal system
215  v_ = getNormal().Cross(u_);
216 }
217 
218 
219 void DetPlane::Print(const Option_t* option) const
220 {
221  printOut<<"DetPlane: "
222  <<"O("<<o_.X()<<", "<<o_.Y()<<", "<<o_.Z()<<") "
223  <<"u("<<u_.X()<<", "<<u_.Y()<<", "<<u_.Z()<<") "
224  <<"v("<<v_.X()<<", "<<v_.Y()<<", "<<v_.Z()<<") "
225  <<"n("<<getNormal().X()<<", "<<getNormal().Y()<<", "<<getNormal().Z()<<") "
226  <<std::endl;
227  if(finitePlane_ != nullptr) {
228  finitePlane_->Print(option);
229  }
230 
231 }
232 
233 
234 /*
235  I could write pages of comments about correct equality checking for
236  floating point numbers, but: When two planes are as close as 10E-5 cm
237  in all nine numbers that define the plane, this will be enough for all
238  practical purposes
239  */
240 bool operator== (const DetPlane& lhs, const DetPlane& rhs){
241  if (&lhs == &rhs)
242  return true;
243  static const double detplaneEpsilon = 1.E-5;
244  if(
245  fabs( (lhs.o_.X()-rhs.o_.X()) ) > detplaneEpsilon ||
246  fabs( (lhs.o_.Y()-rhs.o_.Y()) ) > detplaneEpsilon ||
247  fabs( (lhs.o_.Z()-rhs.o_.Z()) ) > detplaneEpsilon
248  ) return false;
249  else if(
250  fabs( (lhs.u_.X()-rhs.u_.X()) ) > detplaneEpsilon ||
251  fabs( (lhs.u_.Y()-rhs.u_.Y()) ) > detplaneEpsilon ||
252  fabs( (lhs.u_.Z()-rhs.u_.Z()) ) > detplaneEpsilon
253  ) return false;
254  else if(
255  fabs( (lhs.v_.X()-rhs.v_.X()) ) > detplaneEpsilon ||
256  fabs( (lhs.v_.Y()-rhs.v_.Y()) ) > detplaneEpsilon ||
257  fabs( (lhs.v_.Z()-rhs.v_.Z()) ) > detplaneEpsilon
258  ) return false;
259  return true;
260 }
261 
262 bool operator!= (const DetPlane& lhs, const DetPlane& rhs){
263  return !(lhs==rhs);
264 }
265 
266 
267 double DetPlane::distance(const TVector3& point) const {
268  // |(point - o_)*(u_ x v_)|
269  return fabs( (point.X()-o_.X()) * (u_.Y()*v_.Z() - u_.Z()*v_.Y()) +
270  (point.Y()-o_.Y()) * (u_.Z()*v_.X() - u_.X()*v_.Z()) +
271  (point.Z()-o_.Z()) * (u_.X()*v_.Y() - u_.Y()*v_.X()));
272 }
273 
274 double DetPlane::distance(double x, double y, double z) const {
275  // |(point - o_)*(u_ x v_)|
276  return fabs( (x-o_.X()) * (u_.Y()*v_.Z() - u_.Z()*v_.Y()) +
277  (y-o_.Y()) * (u_.Z()*v_.X() - u_.X()*v_.Z()) +
278  (z-o_.Z()) * (u_.X()*v_.Y() - u_.Y()*v_.X()));
279 }
280 
281 
282 TVector2 DetPlane::straightLineToPlane (const TVector3& point, const TVector3& dir) const {
283  TVector3 dirNorm(dir.Unit());
284  TVector3 normal = getNormal();
285  double dirTimesN = dirNorm*normal;
286  if(fabs(dirTimesN)<1.E-6){//straight line is parallel to plane, so return infinity
287  return TVector2(1.E100,1.E100);
288  }
289  double t = 1./dirTimesN * ((o_-point)*normal);
290  return project(point - o_ + t * dirNorm);
291 }
292 
293 
295 void DetPlane::straightLineToPlane(const double& posX, const double& posY, const double& posZ,
296  const double& dirX, const double& dirY, const double& dirZ,
297  double& u, double& v) const {
298 
299  TVector3 W = getNormal();
300  double dirTimesN = dirX*W.X() + dirY*W.Y() + dirZ*W.Z();
301  if(fabs(dirTimesN)<1.E-6){//straight line is parallel to plane, so return infinity
302  u = 1.E100;
303  v = 1.E100;
304  return;
305  }
306  double t = 1./dirTimesN * ((o_.X()-posX)*W.X() +
307  (o_.Y()-posY)*W.Y() +
308  (o_.Z()-posZ)*W.Z());
309 
310  double posOnPlaneX = posX-o_.X() + t*dirX;
311  double posOnPlaneY = posY-o_.Y() + t*dirY;
312  double posOnPlaneZ = posZ-o_.Z() + t*dirZ;
313 
314  u = u_.X()*posOnPlaneX + u_.Y()*posOnPlaneY + u_.Z()*posOnPlaneZ;
315  v = v_.X()*posOnPlaneX + v_.Y()*posOnPlaneY + v_.Z()*posOnPlaneZ;
316 }
317 
318 
319 void DetPlane::rotate(double angle) {
320  TVector3 normal = getNormal();
321  u_.Rotate(angle, normal);
322  v_.Rotate(angle, normal);
323 
324  sane();
325 }
326 
327 
329  o_.SetXYZ(0.,0.,0.);
330  u_.SetXYZ(1.,0.,0.);
331  v_.SetXYZ(0.,1.,0.);
332  finitePlane_.reset();
333 }
334 
335 
336 void DetPlane::Streamer(TBuffer &R__b)
337 {
338  // Stream an object of class genfit::DetPlane.
339  // This is modified from the streamer generated by ROOT in order to
340  // account for the scoped_ptr.
341  //This works around a msvc bug and should be harmless on other platforms
342  typedef ::genfit::DetPlane thisClass;
343  UInt_t R__s, R__c;
344  if (R__b.IsReading()) {
345  Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
346  //TObject::Streamer(R__b);
347  o_.Streamer(R__b);
348  u_.Streamer(R__b);
349  v_.Streamer(R__b);
350  finitePlane_.reset();
351  char flag;
352  R__b >> flag;
353  if (flag) {
354  // Read AbsFinitePlane with correct overload ... see comment
355  // in writing code.
356  TClass* cl = TClass::Load(R__b);
357  AbsFinitePlane *p = (AbsFinitePlane*)(cl->New());
358  cl->Streamer(p, R__b);
359  finitePlane_.reset(p);
360  }
361  R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
362  } else {
363  R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
364  //TObject::Streamer(R__b);
365  o_.Streamer(R__b);
366  u_.Streamer(R__b);
367  v_.Streamer(R__b);
368  if (finitePlane_) {
369  R__b << (char)1;
370 
371  // finitPlane_ is a scoped_ptr, but a shared plane may be
372  // written several times in different places (e.g. it may also
373  // be stored in the measurement class). Since we have no way
374  // of knowing in which other places the same plane is written
375  // (it may even be in another track!), we cannot synchronize
376  // these writes and have to duplicate the SharedPlanePtr's
377  // instead. ROOT caches which pointers were written and read
378  // if one uses 'R__b << p' or equivalent. This can lead to
379  // trouble have no way of synchronizing two reads to
380  // shared_ptr's pointing to the same memory location. But
381  // even if we break the link between the two shared_ptr's,
382  // ROOT will still try to outsmart us, and it will notice that
383  // the finitePlane_ was the same during writing. This will
384  // lead to the same address being attached to two different
385  // scoped_ptrs in reading. Highly undesirable. Since we
386  // cannot know whether other parts of code implicitly or
387  // explicitly rely on TBuffer's caching, we cannot just
388  // disable this caching via R__b.ResetMap() (it's not
389  // documented in any elucidating manner anyway).
390  //
391  // Therefore, we have to write and read the AbsFinitePlane*
392  // manually, bypassing ROOT's caching. In order to do this,
393  // we need the information on the actual type, because
394  // otherwise we couldn't read it back reliably. Finally, the
395  // _working_ means of reading and writing class information
396  // are TClass::Store and TClass::Load, again barely
397  // documented, but if one tries any of the other ways that
398  // suggest themselves, weird breakage will occur.
399  finitePlane_->IsA()->Store(R__b);
400  finitePlane_->Streamer(R__b);
401  } else {
402  R__b << (char)0;
403  }
404  R__b.SetByteCount(R__c, kTRUE);
405  }
406 }
407 } /* End of namespace genfit */