EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LSLTrackRep.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LSLTrackRep.cxx
1 /* Copyright 2008-2009, 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 #include "LSLTrackRep.h"
20 
21 #include <iostream>
22 
23 #include "TMath.h"
24 #include "TMatrixD.h"
25 
26 #include "Nystrom.h"
27 #include "GFAbsRecoHit.h"
28 #include "GFAbsBField.h"
29 #include "LSLEQM.h"
30 #include "AbsNystromEQM.h"
31 
32 LSLTrackRep::LSLTrackRep()
33  : GFAbsTrackRep(5), s(0) ,_acc(1E-2), _adaptive(false)
34 {
35 
36  _eqm=new LSLEQM(0); // this will leak memory in ROOT!
37  fRefPlane=GFDetPlane(TVector3(0,0,0),TVector3(1,0,0),TVector3(0,1,0));
38 
39 }
40 
41 LSLTrackRep::LSLTrackRep(double z, double x, double y,
42  double dxdz, double dydz, double invp,
43  double sigx, double sigy,
44  double sigdxdz, double sigdydz,
45  double siginvp,
46  GFAbsBField* field)
47  : GFAbsTrackRep(5), s(0), _acc(1E-2), _adaptive(false)
48 {
49  s=z;
50  fState[0][0]=x;
51  fState[1][0]=y;
52  fState[2][0]=dxdz;
53  fState[3][0]=dydz;
54  fState[4][0]=invp;
55 
56  fCov[0][0]=sigx;
57  fCov[1][1]=sigy;
58  fCov[2][2]=sigdxdz;
59  fCov[3][3]=sigdydz;
60  fCov[4][4]=siginvp;
61 
62  _eqm=new LSLEQM(field);
63 
64  fRefPlane=GFDetPlane(TVector3(x,y,z),TVector3(1,0,0),TVector3(0,1,0));
65 }
66 
67 LSLTrackRep::LSLTrackRep(const LSLTrackRep& rep)
68  : GFAbsTrackRep(rep)
69 {
70  _eqm=new LSLEQM(0);
71  _acc=rep._acc;
72  _adaptive=rep._adaptive;
73  s=rep.s;
74 }
75 
76 
77 LSLTrackRep::~LSLTrackRep()
78 {
79  if(_eqm!=NULL)delete _eqm;
80 }
81 
82 void
83 LSLTrackRep::init(const TVector3& pos,
84  double dxdz, double dydz, double invp,
85  double sigx, double sigy,
86  double sigdxdz, double sigdydz,
87  double siginvp,
88  GFAbsBField* field)
89 {
90  s=pos.Z();
91  fState[0][0]=pos.X();
92  fState[1][0]=pos.Y();
93  fState[2][0]=dxdz;
94  fState[3][0]=dydz;
95  fState[4][0]=invp;
96 
97  fCov[0][0]=sigx;
98  fCov[1][1]=sigy;
99  fCov[2][2]=sigdxdz;
100  fCov[3][3]=sigdydz;
101  fCov[4][4]=siginvp;
102 
103  if(_eqm!=NULL)delete _eqm;
104  _eqm=new LSLEQM(field);
105 
106  fRefPlane=GFDetPlane(pos,TVector3(1,0,0),TVector3(0,1,0));
107 }
108 
109 
110 void
111 LSLTrackRep::SetBField(GFAbsBField* b)
112 {
113  if(_eqm!=NULL)delete _eqm;
114  _eqm=new LSLEQM(b);
115 }
116 
117 double
118 LSLTrackRep::extrapolate(const GFDetPlane& pl,
119  TMatrixT<double>& statePred)
120 {
121  statePred.ResizeTo(fDimension,1);
122  double sExtrapolateTo=pl.getO().Z();
123  if(sExtrapolateTo<-1000 || sExtrapolateTo>5000)return 0;
124  Nystrom rungeKutta(_eqm);
125  rungeKutta.setAccuracy(_acc);
126  rungeKutta.setAdaptive(_adaptive);
127  //prepare the vectors
128 
129  //std::cout<<"s before extrapolation: "<<s<<std::endl;
130  //std::cout<<"s_to: "<<sExtrapolateTo<<std::endl;
131 
132  TVectorT<double> u(3);u[0]=fState[0][0];u[1]=fState[1][0];u[2]=s;
133  TVectorT<double> uprim(3);uprim[0]=fState[2][0];uprim[1]=fState[3][0];uprim[2]=1.;
134  TVectorT<double> par(1);par[0]=fState[4][0];
135 
136  TVectorT<double> unew(3);unew=u;
137  TVectorT<double> uprimnew(3);uprimnew=uprim;
138  double l=rungeKutta.propagate(s,sExtrapolateTo,
139  u, uprim, par,
140  unew, uprimnew);
141  // write results into statePred
142  statePred[0][0]=unew[0];
143  statePred[1][0]=unew[1];
144  statePred[2][0]=uprimnew[0];
145  statePred[3][0]=uprimnew[1];
146  statePred[4][0]=fState[4][0];
147  //std::cout<<"unew[2]=z="<<unew[2];
148  //std::cout<<"s after extrapolation: "<<s<<std::endl;
149  return l;
150 }
151 
152 /*
153 void LSLTrackRep::extrapolate(const GFDetPlane& pl,
154  const TMatrixT<double>& stateFrom,
155  TMatrixT<double>& stateResult) {
156  stateResult.ResizeTo(5,1);
157  double s=
158  Nystrom rungeKutta(_eqm);
159  //prepare the vectors
160  TVectorT<double> u(3);
161  u[0]=stateFrom[0][0];
162  u[1]=stateFrom[1][0];
163  u[2]=sExtrapolateFrom;
164  TVectorT<double> uprim(3);uprim[0]=stateFrom[2][0];uprim[1]=stateFrom[3][0];uprim[2]=1;
165  TVectorT<double> par(1);par[0]=stateFrom[4][0];
166 
167  TVectorT<double> unew(3);
168  TVectorT<double> uprimnew(3);
169  rungeKutta.propagate(sExtrapolateFrom,sExtrapolateTo,
170  u, uprim, par,
171  unew, uprimnew);
172  // write results into statePred
173  stateResult[0][0]=unew[0];
174  stateResult[1][0]=unew[1];
175  stateResult[2][0]=uprimnew[0];
176  stateResult[3][0]=uprimnew[1];
177  stateResult[4][0]=state[4][0];
178  //std::cout<<"unew[2]=z="<<unew[2];
179 
180 }
181 */
182 
183 double
184 LSLTrackRep::extrapolate(const GFDetPlane& pl,
185  TMatrixT<double>& statePred,
186  TMatrixT<double>& covPred)
187 {
188  statePred.ResizeTo(fDimension,1);
189  covPred.ResizeTo(fDimension,fDimension);
190  TMatrixT<double> jacobian;
191  //std::cout << "Extr from To: " << s << " " << sExtrapolateTo << std::endl;
192  double l=extrapolate(pl,statePred);
193  // covPred=JCovJ^T with J being Jacobian
194  jacobian.ResizeTo(5,5);
195  Jacobian(pl,statePred,jacobian);
196  TMatrixT<double> dummy(fCov,TMatrixT<double>::kMultTranspose,jacobian);
197  covPred=jacobian*dummy;
198  return l;
199 }
200 
201 void
202 LSLTrackRep::extrapolateToPoint(const TVector3& p,
203  TVector3& poca,
204  TVector3& normVec){
205  GFDetPlane plane;
206  int dim = getDim();
207  TMatrixT<double> statePred(dim,1);
208  TMatrixT<double> covPred(dim,dim);
209  plane.setO(p);
210  plane.setU(TVector3(1,0,0));
211  plane.setV(TVector3(0,1,0));
212  extrapolate(plane,statePred,covPred);
213  poca.SetXYZ(statePred[0][0],statePred[1][0],plane.getO().Z());
214  normVec.SetXYZ(statePred[2][0],statePred[3][0],plane.getO().Z());
215 }
216 
217 
218 void
219 LSLTrackRep::stepalong(double h){
220  // create new detplane:
221  GFDetPlane newp(fRefPlane);
222  newp.setO(newp.getO()+TVector3(0,0,h));
223  GFAbsTrackRep::extrapolate(newp);
224 }
225 
226 void
227 LSLTrackRep::Jacobian(const GFDetPlane& pl,
228  const TMatrixT<double>& statePred,
229  TMatrixT<double>& jacResult){
230 
231  TMatrixT<double> difPred(statePred);
232  // do column wise differntiation:
233  for(int icol=0;icol<5;++icol){
234  // choose step
235  double h=TMath::Abs(fState[icol][0])*1.e-4;
236  if(h<1e-13)h=1.e-13;
237  // vary the state
238  fState[icol][0]+=h;
239  extrapolate(pl,difPred);
240  // difference:
241  difPred-=statePred;
242  // remove variation from state
243  fState[icol][0]-=h;
244  // fill jacobian with difference quotient
245  for(int irow=0;irow<5;++irow)jacResult[irow][icol]=difPred[irow][0]/h;
246  }
247 }
248 
249 
250 TVector3
251 LSLTrackRep::getPos(const GFDetPlane& pl)
252 {
253  double z=pl.getO().Z();
254  TMatrixT<double> statePred(fState);
255  GFDetPlane p(TVector3(0,0,z),TVector3(1,0,0),TVector3(0,1,0));
256  extrapolate(p,statePred);
257  return TVector3(statePred[0][0],statePred[1][0],z);
258 }
259 
260 TVector3
261 LSLTrackRep::getMom(const GFDetPlane& pl)
262 {
263  double z=pl.getO().Z();
264  TMatrixT<double> statePred(fState);
265  //statePred.Print();
266  GFDetPlane p(TVector3(0,0,z),TVector3(1,0,0),TVector3(0,1,0));
267  extrapolate(p,statePred);
268  //statePred.Print();
269  TVector3 result(statePred[2][0],statePred[3][0],1);
270  if(TMath::Abs(statePred[4][0])!=0){
271  result.SetMag(1./TMath::Abs(statePred[4][0]));
272  }
273  else result.SetMag(100);
274  if(fInverted)result=(-1.)*result;
275  return result;
276 }
277 
278 void
279 LSLTrackRep::getPosMom(const GFDetPlane& pl,TVector3& pos,TVector3& mom)
280 {
281  double z=pl.getO().Z();
282  TMatrixT<double> statePred(fState);
283  GFDetPlane p(TVector3(0,0,z),TVector3(1,0,0),TVector3(0,1,0));
284  extrapolate(p,statePred);
285  pos.SetXYZ(statePred[0][0],statePred[1][0],z);
286  mom.SetXYZ(statePred[2][0],statePred[3][0],1);
287  if(TMath::Abs(statePred[4][0])!=0){
288  mom.SetMag(1./TMath::Abs(statePred[4][0]));
289  }
290  else mom.SetMag(100);
291  if(fInverted)mom=(-1.)*mom;
292 }
293 
294 
295 TVectorT<double>
296 LSLTrackRep::getGlobal() {// (x,y,z,px,py,pz)
297  TVector3 pos=GFAbsTrackRep::getPos();
298  TVector3 mom=GFAbsTrackRep::getMom();
299  double par[6];
300  par[0]=pos.X();par[1]=pos.Y();par[2]=pos.Z();
301  par[3]=mom.X();par[4]=mom.Y();par[5]=mom.Z();
302  return TVectorT<double>(6,par);
303 }
304 
305 
306 TMatrixT<double>
307 LSLTrackRep::getGlobalCov(){ // covariances
308  TMatrixT<double> L(6,5);
309  double xp=fState[2][0];
310  double yp=fState[3][0];
311  double no=xp*xp+yp*yp+1;
312  double sq=sqrt(no);
313  double sq3inv=1/(sq*sq*sq);
314  double q=getCharge();
315  double p=q/fState[4][0];
316  L[0][0]=1;
317  L[1][1]=1;
318  L[3][2]=p*(1+yp*yp)*sq3inv;
319  L[3][3]=-p*xp*yp*sq3inv;
320  L[3][4]=-p/(2.*fState[4][0])*xp/sq;
321  L[4][2]=p*(1+xp*xp)*sq3inv;
322  L[4][3]=-p*xp*yp*sq3inv;
323  L[4][4]=-p/(2.*fState[4][0])*yp/sq;
324  L[5][2]=-p*xp*sq3inv;
325  L[5][3]=-p*yp*sq3inv;
326  L[5][4]=-p/(2.*fState[4][0])/sq;
327 
328  // calculate new cov;
329  TMatrixT<double> LT(TMatrixD::kTransposed,L);
330  TMatrixT<double> dum(fCov,TMatrixD::kMult,LT);
331  TMatrixT<double> result(L,TMatrixD::kMult,dum);
332 
333  // set sigma_z by hand:
334  result[2][2]=0.01;
335 
336  return result;
337 
338 }
339 
340 ClassImp(LSLTrackRep)