EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GFWireHitPolicy.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GFWireHitPolicy.cxx
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
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 // Please see GFWireHitPolicy.h before using this class.
21 // ----------------------------------------------------------
22 
23 #include "GFWireHitPolicy.h"
24 
25 #include "assert.h"
26 #include <cmath>
27 
28 #include "TMath.h"
29 #include "TVector3.h"
30 
31 #include "GFAbsRecoHit.h"
32 #include "GFException.h"
33 
34 const std::string GFWireHitPolicy::fPolicyName = "GFWireHitPolicy";
35 
36 GFWireHitPolicy::GFWireHitPolicy() :fMaxdistance(1.E50) {;}
37 
38 TMatrixT<double>
40 {
41  TMatrixT<double> returnMat(1,1);
42 
43  checkPlane(hit, plane);
44 
45  // raw x1, y1, z1, x2, y2, z2, rdrift
46  TMatrixT<double> rC = hit->getRawHitCoord();
47 
48  returnMat[0][0] = rC[6][0];
49  return returnMat;
50 }
51 
52 TMatrixT<double>
54 {
55  checkPlane(hit, plane);
56 
57  TMatrixT<double> returnCov(1,1);
58  TMatrixT<double> rawCov = hit->getRawHitCov();
59 
60  returnCov[0][0] = rawCov[6][6];
61 
62  return returnCov;
63 }
64 
65 
66 
68 {
69  // raw x1, y1, z1, x2, y2, z2, rdrift, zreco
70  TMatrixT<double> rC = hit->getRawHitCoord();
71 
72  assert(rC.GetNrows()==7);
73 
74  TVector3 wire1(rC[0][0], rC[1][0], rC[2][0]);
75  TVector3 wire2(rC[3][0], rC[4][0], rC[5][0]);
76  TVector3 wiredirection = wire1 - wire2;
77 
78  TVector3 vaxis = plane.getV();
79  wiredirection.SetMag(1.);
80  vaxis.SetMag(1.);
81 
82  if(fabs(TMath::Abs(wiredirection.Dot(vaxis)) - 1) > 1e-3)
83  {
84 
85  std::cout << "GFWireHitPolicy: plane not valid!!" << std::endl;
86  }
87 }
88 
89 
90 const GFDetPlane&
92 {
93 
94  TMatrixT<double> x=hit->getRawHitCoord();
95  assert(x.GetNrows()==7);
96  TVector3 wire1(x[0][0],x[1][0],x[2][0]);
97  TVector3 wire2(x[3][0],x[4][0],x[5][0]);
98 
99  // distance of one (the first) of the wire extremities from the plane
100  Double_t d_from_refplane = fDetPlane.dist(wire1).Mag();
101  if(d_from_refplane < 1e-5) return fDetPlane;
102 
103 
104  // point of closest approach
105  TVector3 poca, poca_onwire, dirInPoca;
106 
107  rep->extrapolateToLine(wire1, wire2, poca, dirInPoca, poca_onwire);
108 
109 
110  Double_t distance;
111  distance = TMath::Sqrt(fabs(((wire1-poca).Mag2()*(wire2-wire1).Mag2()-pow((wire1-poca).Dot(wire2-wire1),2))/(wire2-wire1).Mag2()));
112 
113  // check poca inside tube
114  if(distance > fMaxdistance) {
115  GFException exc("distance poca-wire > maxdistance", __LINE__,__FILE__);
116  throw exc;
117  }
118 
119  // find plane
120  // unitary vector along distance
121  // poca (on track), poca_onwire (on wire)
122  TVector3 fromwiretoextr = poca - poca_onwire;
123  fromwiretoextr.SetMag(1.);
124  // unitary vector along the wire
125  TVector3 wiredirection = wire2 - wire1;
126  wiredirection.SetMag(1.);
127 
128  // check orthogonality
129  if(fabs(fromwiretoextr * wiredirection) > 1e-3) {
130  GFException exc("fromwiretoextr*wiredirection > 1e-3", __LINE__,__FILE__);
131  throw exc;
132  }
133 
134  TVector3 U;
135  U = fromwiretoextr;
136  TVector3 V;
137  V = wiredirection;
138  U.SetMag(1.);
139  V.SetMag(1.);
140 
141  TVector3 O = poca_onwire;
142 
143  fDetPlane = GFDetPlane(O, U, V);
144 
145  return fDetPlane;
146 }
147