FairRoot/PandaRoot
RKTrackRep.cxx
Go to the documentation of this file.
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 
20 /* The Runge Kutta implementation stems from GEANT3 originally (R. Brun et al.)
21  Porting to C goes back to Igor Gavrilenko @ CERN
22  The code was taken from the Phast analysis package of the COMPASS experiment
23  (Sergei Gerrassimov @ CERN)
24 */
25 
26 #include"RKTrackRep.h"
27 #include <iostream>
28 #include <sstream>
29 #include <iomanip>
30 #include"assert.h"
31 #include "stdlib.h"
32 #include"math.h"
33 #include"TMath.h"
34 #include"TGeoManager.h"
35 #include"TDatabasePDG.h"
36 #include"GFException.h"
37 #include"GFFieldManager.h"
38 #include"GFMaterialEffects.h"
39 
40 #define MINSTEP 0.001 // minimum step [cm] for Runge Kutta and iteration to POCA
41 
42 
43 void RKTrackRep::setData(const TMatrixT<double>& st, const GFDetPlane& pl, const TMatrixT<double>* cov, const TMatrixT<double>* aux){
44  if(aux != NULL) {
45  fCacheSpu = (*aux)(0,0);
46  } else {
47  if(pl!=fCachePlane){
48  std::cerr << "RKTrackRep::setData() - a fatal error ocurred! It was called with a reference plane which is not the same as the one from the last extrapolate(plane,state,cov)-> abort in line " << __LINE__ << std::endl;
49  throw;
50  }
51  }
52  GFAbsTrackRep::setData(st,pl,cov);
53  if (fCharge*fState[0][0] < 0) fCharge *= -1; // set charge accordingly! (fState[0][0] = q/p)
54  fSpu=fCacheSpu;
55 }
56 
57 const TMatrixT<double>* RKTrackRep::getAuxInfo(const GFDetPlane& pl) {
58 
59  if(pl!=fCachePlane) {
60  std::cerr << "RKTrackRep::getAuxInfo() - Fatal error: Trying to get auxillary information with planes mismatch (Information returned does not belong to requested plane)! -> abort in line " << __LINE__ << std::endl;
61  throw;
62  }
63  fAuxInfo.ResizeTo(1,1);
64  fAuxInfo(0,0) = fCacheSpu;
65  return &fAuxInfo;
66 
67 }
68 
69 RKTrackRep::~RKTrackRep(){
70  //GFMaterialEffects is now a Singleton
71 }
72 
73 
74 RKTrackRep::RKTrackRep() : GFAbsTrackRep(5), fDirection(true), fPdg(0), fMass(0.), fCharge(-1), fCachePlane(), fCacheSpu(1), fSpu(1), fAuxInfo(1,1) {
75 }
76 
77 
78 RKTrackRep::RKTrackRep(const TVector3& pos,
79  const TVector3& mom,
80  const TVector3& poserr,
81  const TVector3& momerr,
82  const int& PDGCode) :
83  GFAbsTrackRep(5), fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1) {
84  setPDG(PDGCode); // also sets charge and mass
85 
86 
87  fRefPlane.setO(pos);
88  fRefPlane.setNormal(mom);
89 
90  fState[0][0]=fCharge/mom.Mag();
91  //u' and v'
92  fState[1][0]=0.0;
93  fState[2][0]=0.0;
94  //u and v
95  fState[3][0]=0.0;
96  fState[4][0]=0.0;
97  //spu
98  fSpu=1.;
99 
100  TVector3 o=fRefPlane.getO();
101  TVector3 u=fRefPlane.getU();
102  TVector3 v=fRefPlane.getV();
103  TVector3 w=u.Cross(v);
104  double pu=0.;
105  double pv=0.;
106  double pw=mom.Mag();
107 
108  fCov[3][3] = poserr.X()*poserr.X() * u.X()*u.X() +
109  poserr.Y()*poserr.Y() * u.Y()*u.Y() +
110  poserr.Z()*poserr.Z() * u.Z()*u.Z();
111  fCov[4][4] = poserr.X()*poserr.X() * v.X()*v.X() +
112  poserr.Y()*poserr.Y() * v.Y()*v.Y() +
113  poserr.Z()*poserr.Z() * v.Z()*v.Z();
114  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
115  (mom.X()*mom.X() * momerr.X()*momerr.X()+
116  mom.Y()*mom.Y() * momerr.Y()*momerr.Y()+
117  mom.Z()*mom.Z() * momerr.Z()*momerr.Z());
118  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*momerr.X()*momerr.X() +
119  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
120  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*momerr.Z()*momerr.Z();
121  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*momerr.X()*momerr.X() +
122  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
123  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*momerr.Z()*momerr.Z();
124 
125 }
126 
127 
128 RKTrackRep::RKTrackRep(const GFTrackCand* aGFTrackCandPtr) :
129  GFAbsTrackRep(5), fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1) {
130  setPDG(aGFTrackCandPtr->getPdgCode()); // also sets charge and mass
131 
132  TVector3 mom = aGFTrackCandPtr->getDirSeed();
133  fRefPlane.setO(aGFTrackCandPtr->getPosSeed());
134  fRefPlane.setNormal(mom);
135  double pw=mom.Mag();
136  fState[0][0]=fCharge/pw;
137  //u' and v'
138  fState[1][0]=0.0;
139  fState[2][0]=0.0;
140  //u and v
141  fState[3][0]=0.0;
142  fState[4][0]=0.0;
143  //spu
144  fSpu=1.;
145 
146  TVector3 o=fRefPlane.getO();
147  TVector3 u=fRefPlane.getU();
148  TVector3 v=fRefPlane.getV();
149  TVector3 w=u.Cross(v);
150  double pu=0.;
151  double pv=0.;
152 
153  TVector3 poserr = aGFTrackCandPtr->getPosError();
154  TVector3 momerr = aGFTrackCandPtr->getDirError();
155  fCov[3][3] = poserr.X()*poserr.X() * u.X()*u.X() +
156  poserr.Y()*poserr.Y() * u.Y()*u.Y() +
157  poserr.Z()*poserr.Z() * u.Z()*u.Z();
158  fCov[4][4] = poserr.X()*poserr.X() * v.X()*v.X() +
159  poserr.Y()*poserr.Y() * v.Y()*v.Y() +
160  poserr.Z()*poserr.Z() * v.Z()*v.Z();
161  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
162  (mom.X()*mom.X() * momerr.X()*momerr.X()+
163  mom.Y()*mom.Y() * momerr.Y()*momerr.Y()+
164  mom.Z()*mom.Z() * momerr.Z()*momerr.Z());
165  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*momerr.X()*momerr.X() +
166  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
167  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*momerr.Z()*momerr.Z();
168  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*momerr.X()*momerr.X() +
169  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*momerr.Y()*momerr.Y() +
170  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*momerr.Z()*momerr.Z();
171 
172 }
173 
174 RKTrackRep::RKTrackRep(const TVector3& pos,
175  const TVector3& mom,
176  const int& PDGCode) :
177  GFAbsTrackRep(5),fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1){
178  setPDG(PDGCode); // also sets charge and mass
179 
180 
181  fRefPlane.setO(pos);
182  fRefPlane.setNormal(mom);
183 
184  fState[0][0]=fCharge/mom.Mag();
185  //u' and v'
186  fState[1][0]=0.0;
187  fState[2][0]=0.0;
188  //u and v
189  fState[3][0]=0.0;
190  fState[4][0]=0.0;
191  //spu
192  fSpu=1.;
193 
194  TVector3 o=fRefPlane.getO();
195  TVector3 u=fRefPlane.getU();
196  TVector3 v=fRefPlane.getV();
197  TVector3 w=u.Cross(v);
198  double pu=0.;
199  double pv=0.;
200  double pw=mom.Mag();
201 
202  static const TVector3 stdPosErr(1.,1.,1.);
203  static const TVector3 stdMomErr(1.,1.,1.);
204 
205  fCov[3][3] = stdPosErr.X()*stdPosErr.X() * u.X()*u.X() +
206  stdPosErr.Y()*stdPosErr.Y() * u.Y()*u.Y() +
207  stdPosErr.Z()*stdPosErr.Z() * u.Z()*u.Z();
208  fCov[4][4] = stdPosErr.X()*stdPosErr.X() * v.X()*v.X() +
209  stdPosErr.Y()*stdPosErr.Y() * v.Y()*v.Y() +
210  stdPosErr.Z()*stdPosErr.Z() * v.Z()*v.Z();
211  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
212  (mom.X()*mom.X() * stdMomErr.X()*stdMomErr.X()+
213  mom.Y()*mom.Y() * stdMomErr.Y()*stdMomErr.Y()+
214  mom.Z()*mom.Z() * stdMomErr.Z()*stdMomErr.Z());
215  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*stdMomErr.X()*stdMomErr.X() +
216  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*stdMomErr.Y()*stdMomErr.Y() +
217  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*stdMomErr.Z()*stdMomErr.Z();
218  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*stdMomErr.X()*stdMomErr.X() +
219  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*stdMomErr.Y()*stdMomErr.Y() +
220  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*stdMomErr.Z()*stdMomErr.Z();
221 
222 }
223 
224 RKTrackRep::RKTrackRep(const GFDetPlane& pl,
225  const TVector3& mom,
226  const int& PDGCode) :
227  GFAbsTrackRep(5),fDirection(true), fCachePlane(), fCacheSpu(1), fAuxInfo(1,1){
228  setPDG(PDGCode); // also sets charge and mass
229 
230 
231  fRefPlane = pl;
232  TVector3 o=fRefPlane.getO();
233  TVector3 u=fRefPlane.getU();
234  TVector3 v=fRefPlane.getV();
235  TVector3 w=u.Cross(v);
236 
237  double pu=mom*u;
238  double pv=mom*v;
239  double pw=mom*w;
240 
241  fState[0][0]=fCharge/mom.Mag();
242  //u' and v'
243  fState[1][0]=pu/pw;
244  fState[2][0]=pv/pw;
245  //u and v
246  fState[3][0]=0.0;
247  fState[4][0]=0.0;
248  //spu
249  fSpu=1.;
250 
251  static const TVector3 stdPosErr2(1.,1.,1.);
252  static const TVector3 stdMomErr2(1.,1.,1.);
253 
254  fCov[3][3] = stdPosErr2.X()*stdPosErr2.X() * u.X()*u.X() +
255  stdPosErr2.Y()*stdPosErr2.Y() * u.Y()*u.Y() +
256  stdPosErr2.Z()*stdPosErr2.Z() * u.Z()*u.Z();
257  fCov[4][4] = stdPosErr2.X()*stdPosErr2.X() * v.X()*v.X() +
258  stdPosErr2.Y()*stdPosErr2.Y() * v.Y()*v.Y() +
259  stdPosErr2.Z()*stdPosErr2.Z() * v.Z()*v.Z();
260  fCov[0][0] = fCharge*fCharge/pow(mom.Mag(),6.) *
261  (mom.X()*mom.X() * stdMomErr2.X()*stdMomErr2.X()+
262  mom.Y()*mom.Y() * stdMomErr2.Y()*stdMomErr2.Y()+
263  mom.Z()*mom.Z() * stdMomErr2.Z()*stdMomErr2.Z());
264  fCov[1][1] = pow((u.X()/pw - w.X()*pu/(pw*pw)),2.)*stdMomErr2.X()*stdMomErr2.X() +
265  pow((u.Y()/pw - w.Y()*pu/(pw*pw)),2.)*stdMomErr2.Y()*stdMomErr2.Y() +
266  pow((u.Z()/pw - w.Z()*pu/(pw*pw)),2.)*stdMomErr2.Z()*stdMomErr2.Z();
267  fCov[2][2] = pow((v.X()/pw - w.X()*pv/(pw*pw)),2.)*stdMomErr2.X()*stdMomErr2.X() +
268  pow((v.Y()/pw - w.Y()*pv/(pw*pw)),2.)*stdMomErr2.Y()*stdMomErr2.Y() +
269  pow((v.Z()/pw - w.Z()*pv/(pw*pw)),2.)*stdMomErr2.Z()*stdMomErr2.Z();
270 
271 }
272 
273 
274 void RKTrackRep::setPDG(int i){
275  fPdg = i;
276  TParticlePDG * part = TDatabasePDG::Instance()->GetParticle(fPdg);
277  if(part == 0){
278  std::cerr << "RKTrackRep::setPDG particle " << i
279  << " not known to TDatabasePDG -> abort" << std::endl;
280  exit(1);
281  }
282  fMass = part->Mass();
283  fCharge = part->Charge()/(3.);
284 }
285 
286 
287 TVector3 RKTrackRep::getPos(const GFDetPlane& pl){
288  if(pl!=fRefPlane){
289  TMatrixT<double> s(5,1);
290  extrapolate(pl,s);
291  return pl.getO()+s[3][0]*pl.getU()+s[4][0]*pl.getV();
292  }
293  return fRefPlane.getO()+fState[3][0]*fRefPlane.getU()+fState[4][0]*fRefPlane.getV();
294 }
295 
296 
297 TVector3 RKTrackRep::getMom(const GFDetPlane& pl){
298  TMatrixT<double> statePred(fState);
299  TVector3 retmom;
300  if(pl!=fRefPlane) {
301  extrapolate(pl,statePred);
302  retmom = fCacheSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
303  }
304  else{
305  retmom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
306  }
307  retmom.SetMag(1./fabs(statePred[0][0]));
308  return retmom;
309 }
310 
311 
312 void RKTrackRep::getPosMom(const GFDetPlane& pl,TVector3& pos,
313  TVector3& mom){
314  TMatrixT<double> statePred(fState);
315  if(pl!=fRefPlane) {
316  extrapolate(pl,statePred);
317  mom = fCacheSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
318  }
319  else{
320  mom = fSpu*(pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV());
321  }
322  mom.SetMag(1./fabs(statePred[0][0]));
323  pos = pl.getO()+(statePred[3][0]*pl.getU())+(statePred[4][0]*pl.getV());
324 }
325 
326 
327 
328 
329 void RKTrackRep::extrapolateToPoint(const TVector3& pos,
330  TVector3& poca,
331  TVector3& dirInPoca){
332 
333  static const int maxIt(30);
334 
335  TVector3 o=fRefPlane.getO();
336  TVector3 u=fRefPlane.getU();
337  TVector3 v=fRefPlane.getV();
338  TVector3 w=u.Cross(v);
339 
340  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
341  pTilde.SetMag(1.);
342 
343  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
344 
345  TMatrixT<double> state7(7,1);
346  state7[0][0] = point.X();
347  state7[1][0] = point.Y();
348  state7[2][0] = point.Z();
349  state7[3][0] = pTilde.X();
350  state7[4][0] = pTilde.Y();
351  state7[5][0] = pTilde.Z();
352  state7[6][0] = fState[0][0];
353 
354  double coveredDistance(0.);
355 
356  GFDetPlane pl;
357  int iterations(0);
358 
359  while(true){
360  pl.setON(pos,TVector3(state7[3][0],state7[4][0],state7[5][0]));
361  coveredDistance = this->Extrap(pl,&state7);
362 
363  if(fabs(coveredDistance)<MINSTEP) break;
364  if(++iterations == maxIt) {
365  GFException exc("RKTrackRep::extrapolateToPoint ==> extrapolation to point failed, maximum number of iterations reached",__LINE__,__FILE__);
366  throw exc;
367  }
368  }
369  poca.SetXYZ(state7[0][0],state7[1][0],state7[2][0]);
370  dirInPoca.SetXYZ(state7[3][0],state7[4][0],state7[5][0]);
371 }
372 
373 
374 
375 
376 TVector3 RKTrackRep::poca2Line(const TVector3& extr1,const TVector3& extr2,const TVector3& point) const {
377 
378  TVector3 theWire = extr2-extr1;
379  if(theWire.Mag()<1.E-8){
380  GFException exc("RKTrackRep::poca2Line ==> try to find poca between line and point, but the line is really just a point",__LINE__,__FILE__);
381  throw exc;
382  }
383  double t = 1./(theWire*theWire)*(point*theWire+extr1*extr1-extr1*extr2);
384  return (extr1+t*theWire);
385 }
386 
387 
388 
389 
390 void RKTrackRep::extrapolateToLine(const TVector3& point1,
391  const TVector3& point2,
392  TVector3& poca,
393  TVector3& dirInPoca,
394  TVector3& poca_onwire){
395  static const int maxIt(30);
396 
397  TVector3 o=fRefPlane.getO();
398  TVector3 u=fRefPlane.getU();
399  TVector3 v=fRefPlane.getV();
400  TVector3 w=u.Cross(v);
401 
402  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
403  pTilde.SetMag(1.);
404 
405  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
406 
407  TMatrixT<double> state7(7,1);
408  state7[0][0] = point.X();
409  state7[1][0] = point.Y();
410  state7[2][0] = point.Z();
411  state7[3][0] = pTilde.X();
412  state7[4][0] = pTilde.Y();
413  state7[5][0] = pTilde.Z();
414  state7[6][0] = fState[0][0];
415 
416  double coveredDistance(0.);
417 
418  GFDetPlane pl;
419  int iterations(0);
420 
421  while(true){
422  pl.setO(point1);
423  TVector3 currentDir(state7[3][0],state7[4][0],state7[5][0]);
424  pl.setU(currentDir.Cross(point2-point1));
425  pl.setV(point2-point1);
426  coveredDistance = this->Extrap(pl,&state7);
427 
428  if(fabs(coveredDistance)<MINSTEP) break;
429  if(++iterations == maxIt) {
430  GFException exc("RKTrackRep::extrapolateToLine ==> extrapolation to line failed, maximum number of iterations reached",__LINE__,__FILE__);
431  throw exc;
432  }
433  }
434  poca.SetXYZ(state7[0][0],state7[1][0],state7[2][0]);
435  dirInPoca.SetXYZ(state7[3][0],state7[4][0],state7[5][0]);
436  poca_onwire = poca2Line(point1,point2,poca);
437 }
438 
439 
440 
441 
442 double RKTrackRep::extrapolate(const GFDetPlane& pl,
443  TMatrixT<double>& statePred,
444  TMatrixT<double>& covPred){
445 
446  TMatrixT<double> cov7x7(7,7);
447  TMatrixT<double> J_pM(7,5);
448 
449  TVector3 o=fRefPlane.getO();
450  TVector3 u=fRefPlane.getU();
451  TVector3 v=fRefPlane.getV();
452  TVector3 w=u.Cross(v);
453 
454  J_pM[0][3] = u.X();J_pM[0][4]=v.X(); // dx/du
455  J_pM[1][3] = u.Y();J_pM[1][4]=v.Y();
456  J_pM[2][3] = u.Z();J_pM[2][4]=v.Z();
457 
458  TVector3 pTilde = fSpu * (w + fState[1][0] * u + fState[2][0] * v);
459  double pTildeMag = pTilde.Mag();
460 
461  //J_pM matrix is d(x,y,z,ax,ay,az,q/p) / d(q/p,u',v',u,v)
462 
463  // da_x/du'
464  J_pM[3][1] = fSpu/pTildeMag*(u.X()-pTilde.X()/(pTildeMag*pTildeMag)*u*pTilde);
465  J_pM[4][1] = fSpu/pTildeMag*(u.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*u*pTilde);
466  J_pM[5][1] = fSpu/pTildeMag*(u.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*u*pTilde);
467  // da_x/dv'
468  J_pM[3][2] = fSpu/pTildeMag*(v.X()-pTilde.X()/(pTildeMag*pTildeMag)*v*pTilde);
469  J_pM[4][2] = fSpu/pTildeMag*(v.Y()-pTilde.Y()/(pTildeMag*pTildeMag)*v*pTilde);
470  J_pM[5][2] = fSpu/pTildeMag*(v.Z()-pTilde.Z()/(pTildeMag*pTildeMag)*v*pTilde);
471  // dqOp/dqOp
472  J_pM[6][0] = 1.;
473 
474  TMatrixT<double> J_pM_transp(J_pM);
475  J_pM_transp.T();
476 
477  cov7x7 = J_pM*(fCov*J_pM_transp);
478 
479  TVector3 pos = o + fState[3][0]*u + fState[4][0]*v;
480  TMatrixT<double> state7(7,1);
481  state7[0][0] = pos.X();
482  state7[1][0] = pos.Y();
483  state7[2][0] = pos.Z();
484  state7[3][0] = pTilde.X()/pTildeMag;;
485  state7[4][0] = pTilde.Y()/pTildeMag;;
486  state7[5][0] = pTilde.Z()/pTildeMag;;
487  state7[6][0] = fState[0][0];
488 
489  double coveredDistance = this->Extrap(pl,&state7,&cov7x7);
490 
491  TVector3 O = pl.getO();
492  TVector3 U = pl.getU();
493  TVector3 V = pl.getV();
494  TVector3 W = pl.getNormal();
495 
496  double X = state7[0][0];
497  double Y = state7[1][0];
498  double Z = state7[2][0];
499  double AX = state7[3][0];
500  double AY = state7[4][0];
501  double AZ = state7[5][0];
502  double QOP = state7[6][0];
503  TVector3 A(AX,AY,AZ);
504  TVector3 Point(X,Y,Z);
505  TMatrixT<double> J_Mp(5,7);
506 
507  // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,ax,ay,az,q/p)
508  J_Mp[0][6] = 1.;
509  //du'/da_x
510  double AtW = A*W;
511  J_Mp[1][3] = (U.X()*(AtW)-W.X()*(A*U))/(AtW*AtW);
512  J_Mp[1][4] = (U.Y()*(AtW)-W.Y()*(A*U))/(AtW*AtW);
513  J_Mp[1][5] = (U.Z()*(AtW)-W.Z()*(A*U))/(AtW*AtW);
514  //dv'/da_x
515  J_Mp[2][3] = (V.X()*(AtW)-W.X()*(A*V))/(AtW*AtW);
516  J_Mp[2][4] = (V.Y()*(AtW)-W.Y()*(A*V))/(AtW*AtW);
517  J_Mp[2][5] = (V.Z()*(AtW)-W.Z()*(A*V))/(AtW*AtW);
518  //du/dx
519  J_Mp[3][0] = U.X();
520  J_Mp[3][1] = U.Y();
521  J_Mp[3][2] = U.Z();
522  //dv/dx
523  J_Mp[4][0] = V.X();
524  J_Mp[4][1] = V.Y();
525  J_Mp[4][2] = V.Z();
526 
527  TMatrixT<double> J_Mp_transp(J_Mp);
528  J_Mp_transp.T();
529 
530  covPred.ResizeTo(5,5);
531 
532  covPred = J_Mp*(cov7x7*J_Mp_transp);
533 
534  statePred.ResizeTo(5,1);
535  statePred[0][0] = QOP;
536  statePred[1][0] = (A*U)/(A*W);
537  statePred[2][0] = (A*V)/(A*W);
538  statePred[3][0] = (Point-O)*U;
539  statePred[4][0] = (Point-O)*V;
540 
541  fCachePlane = pl;
542  fCacheSpu = (A*W)/fabs(A*W);
543 
544  return coveredDistance;
545 }
546 
547 
548 
549 
550 double RKTrackRep::extrapolate(const GFDetPlane& pl,
551  TMatrixT<double>& statePred){
552 
553  TVector3 o=fRefPlane.getO();
554  TVector3 u=fRefPlane.getU();
555  TVector3 v=fRefPlane.getV();
556  TVector3 w=u.Cross(v);
557 
558 
559  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
560  double pTildeMag = pTilde.Mag();
561 
562  TVector3 pos = o + fState[3][0]*u + fState[4][0]*v;
563 
564  TMatrixT<double> state7(7,1);
565  state7[0][0] = pos.X();
566  state7[1][0] = pos.Y();
567  state7[2][0] = pos.Z();
568  state7[3][0] = pTilde.X()/pTildeMag;
569  state7[4][0] = pTilde.Y()/pTildeMag;
570  state7[5][0] = pTilde.Z()/pTildeMag;
571  state7[6][0] = fState[0][0];
572 
573  TVector3 O = pl.getO();
574  TVector3 U = pl.getU();
575  TVector3 V = pl.getV();
576  TVector3 W = pl.getNormal();
577 
578  double coveredDistance = this->Extrap(pl,&state7);
579 
580  double X = state7[0][0];
581  double Y = state7[1][0];
582  double Z = state7[2][0];
583  double AX = state7[3][0];
584  double AY = state7[4][0];
585  double AZ = state7[5][0];
586  double QOP = state7[6][0];
587  TVector3 A(AX,AY,AZ);
588  TVector3 Point(X,Y,Z);
589 
590  statePred.ResizeTo(5,1);
591  statePred[0][0] = QOP;
592  statePred[1][0] = (A*U)/(A*W);
593  statePred[2][0] = (A*V)/(A*W);
594  statePred[3][0] = (Point-O)*U;
595  statePred[4][0] = (Point-O)*V;
596 
597  return coveredDistance;
598 }
599 
600 
601 
602 
603 //
604 // Runge-Kutta method for tracking a particles through a magnetic field.
605 // Uses Nystroem algorithm (See Handbook Nat. Bur. of Standards, procedure 25.5.20)
606 //
607 // Input parameters:
608 // SU - plane parameters
609 // SU[0] - direction cosines normal to surface Ex
610 // SU[1] - ------- Ey
611 // SU[2] - ------- Ez; Ex*Ex+Ey*Ey+Ez*Ez=1
612 // SU[3] - distance to surface from (0,0,0) > 0 cm
613 //
614 // ND - number of variables for derivatives calculation
615 // P - initial parameters (coordinates(cm), direction cosines,
616 // charge/momentum (Gev-1) and derivatives this parameters (8x7)
617 //
618 // X Y Z Ax Ay Az q/P
619 // P[ 0] P[ 1] P[ 2] P[ 3] P[ 4] P[ 5] P[ 6]
620 //
621 // dX/dp dY/dp dZ/dp dAx/dp dAy/dp dAz/dp d(q/P)/dp*P[6]
622 // P[ 7] P[ 8] P[ 9] P[10] P[11] P[12] P[13] d()/dp1
623 //
624 // P[14] P[15] P[16] P[17] P[18] P[19] P[20] d()/dp2
625 // ............................................................................ d()/dpND
626 //
627 // Output parameters:
628 //
629 // P - output parameters and derivatives after propagation in magnetic field
630 // defined by Mfield (KGauss)
631 // Where a Mfield(R,H) - is interface to magnetic field information
632 // input R[ 0],R[ 1],R[ 2] - X , Y and Z of the track
633 // output H[ 0],H[ 1],H[ 2] - Hx , Hy and Hz of the magnetic field
634 // H[ 3],H[ 4],H[ 5] - dHx/dx, dHx/dy and dHx/dz //
635 // H[ 6],H[ 7],H[ 8] - dHy/dx, dHy/dy and dHy/dz // (not used)
636 // H[ 9],H[10],H[11] - dHz/dx, dHz/dy and dHz/dz //
637 //
638 // Authors: R.Brun, M.Hansroul, V.Perevoztchikov (Geant3)
639 //
640 bool RKTrackRep::RKutta (const GFDetPlane& plane,
641  double* P,
642  double& coveredDistance,
643  std::vector<TVector3>& points,
644  std::vector<double>& pointPaths,
645  const double& , // currently not used // maxLen /
646  bool calcCov) const {
647 
648  static const double EC = .000149896229; // c/(2*10^12) resp. c/2Tera
649  static const double DLT = .0002; // max. deviation for approximation-quality test
650  static const double DLT32 = DLT/32.; //
651  static const double P3 = 1./3.; // 1/3
652  static const double Smax = 100.; // max. step allowed > 0
653  static const double Wmax = 3000.; // max. way allowed
654  static const double Pmin = 4.E-3; // minimum momentum for propagation [GeV]
655  static const int ND = 56; // number of variables for derivatives calculation
656  static const int ND1 = ND-7; // = 49
657  double* R = &P[0]; // Start coordinates in cm ( x, y, z)
658  double* A = &P[3]; // Start directions (ax, ay, az); ax^2+ay^2+az^2=1
659  double SA[3] = {0.,0.,0.}; // Start directions derivatives
660  double Pinv = P[6]*EC; // P[6] is charge/momentum in e/(Gev/c)
661  double Way = 0.; // Total way of the trajectory
662  double Way2 = 0.; // Total way of the trajectory with correct signs
663  bool error = false; // Error of propogation
664  bool stopBecauseOfMaterial = false; // does not go through main loop again when stepsize is reduced by stepper
665 
666  points.clear();
667  pointPaths.clear();
668  //std::cout<<"coords "<<R[0]<<" "<<R[1]<<" "<<R[2]<< std::endl;
669  //std::cout<<"R "<<sqrt(pow(R[0],2)+pow(R[1],2))<< std::endl;
670  //std::cout<<"momentum "<<fabs(fCharge/P[6])<< std::endl;
671  if(fabs(fCharge/P[6])<Pmin){
672  std::cerr << "RKTrackRep::RKutta ==> momentum too low: " << fabs(fCharge/P[6])*1000. << " MeV" << std::endl;
673  return (false);
674  }
675 
676  double SU[4];
677  TVector3 O = plane.getO();
678  TVector3 W = plane.getNormal();
679  if(W*O > 0){ // make SU vector point away from origin
680  SU[0] = W.X();
681  SU[1] = W.Y();
682  SU[2] = W.Z();
683  }
684  else{
685  SU[0] = -1.*W.X();
686  SU[1] = -1.*W.Y();
687  SU[2] = -1.*W.Z();
688  }
689  SU[3] = plane.distance(0.,0.,0.);
690 
691  //
692  // Step estimation until surface
693  //
694  double Step,An,Dist,Dis,S,Sl=0;
695 
696  points.push_back(TVector3(R[0],R[1],R[2]));
697  pointPaths.push_back(0.);
698 
699  An=A[0]*SU[0]+A[1]*SU[1]+A[2]*SU[2]; // An = A * N; component of A normal to surface
700  if(fabs(An) < 1.E-6) {
701  std::cerr << "RKTrackRep::RKutta ==> cannot propagate perpendicular to plane " << std::endl;
702  return false; // no propagation if A perpendicular to surface
703  }
704  if( plane.inActive(TVector3(R[0],R[1],R[2]),TVector3(A[0],A[1],A[2]))) { // if direction is not pointing to active part of surface
705  Dist=SU[3]-R[0]*SU[0]-R[1]*SU[1]-R[2]*SU[2]; // Distance between start coordinates and surface
706  Step=Dist/An;
707  }
708  else{ // make sure to extrapolate towards surface
709  if( (O.X()-R[0])*A[0] + (O.Y()-R[1])*A[1] + (O.Z()-R[2])*A[2] >0 ){ // if direction A pointing from start coordinates R towards surface
710  Dist = sqrt((R[0]-O.X())*(R[0]-O.X())+ // |R-O|; Distance between start coordinates and origin of surface
711  (R[1]-O.Y())*(R[1]-O.Y())+
712  (R[2]-O.Z())*(R[2]-O.Z()));
713  }
714  else{ // if direction pointing away from surface
715  Dist = -1.*sqrt((R[0]-O.X())*(R[0]-O.X())+
716  (R[1]-O.Y())*(R[1]-O.Y())+
717  (R[2]-O.Z())*(R[2]-O.Z()));
718  }
719  Step=Dist;
720  }
721 
722  if(fabs(Step)>Wmax) {
723  std::cerr<<"RKTrackRep::RKutta ==> Too long extrapolation requested : "<<Step<<" cm !"<<std::endl;
724  std::cerr<<"X = "<<R[0]<<" Y = "<<R[1]<<" Z = "<<R[2]
725  <<" COSx = "<<A[0]<<" COSy = "<<A[1]<<" COSz = "<<A[2]<<std::endl;
726  std::cout<<"Destination X = "<<SU[0]*SU[3]<<std::endl;
727  return(false);
728  }
729 
730  // reduce maximum stepsize S to Smax
731  Step>Smax ? S=Smax : Step<-Smax ? S=-Smax : S=Step;
732 
733  //
734  // Main cycle of Runge-Kutta method
735  //
736 
737  //for saving points only when the direction didnt change
738  int Ssign=1;
739  if(S<0) Ssign = -1;
740 
741  while(fabs(Step)>MINSTEP && !stopBecauseOfMaterial) {
742 
743  // call stepper and reduce stepsize
744  double stepperLen;
745  stepperLen = GFMaterialEffects::getInstance()->stepper(fabs(S),
746  R[0],R[1],R[2],
747  Ssign*A[0],Ssign*A[1],Ssign*A[2],
748  fabs(fCharge/P[6]),
749  fPdg);
750  if (stepperLen<MINSTEP) stepperLen=MINSTEP; // prevents tiny stepsizes that can slow down the program
751  if (S > stepperLen) {
752  S = stepperLen;
753  stopBecauseOfMaterial = true;
754  }
755  else if (S < -stepperLen) {
756  S = -stepperLen;
757  stopBecauseOfMaterial = true;
758  }
759 
760  double H0[12],H1[12],H2[12],r[3];
761  double S3=P3*S, S4=.25*S, PS2=Pinv*S;
762 
763  //
764  // First point
765  //
766  r[0]=R[0] ; r[1]=R[1] ; r[2]=R[2] ;
767  TVector3 pos(r[0],r[1],r[2]); // vector of start coordinates R0 (x, y, z)
768  TVector3 H0vect = GFFieldManager::getFieldVal(pos); // magnetic field in 10^-4 T = kGauss
769  H0[0]=PS2*H0vect.X(); H0[1]=PS2*H0vect.Y(); H0[2]=PS2*H0vect.Z(); // H0 is PS2*(Hx, Hy, Hz) @ R0
770  double A0=A[1]*H0[2]-A[2]*H0[1], B0=A[2]*H0[0]-A[0]*H0[2], C0=A[0]*H0[1]-A[1]*H0[0]; // (ax, ay, az) x H0
771  double A2=A[0]+A0 , B2=A[1]+B0 , C2=A[2]+C0 ; // (A0, B0, C0) + (ax, ay, az)
772  double A1=A2+A[0] , B1=B2+A[1] , C1=C2+A[2] ; // (A0, B0, C0) + 2*(ax, ay, az)
773 
774  //
775  // Second point
776  //
777  r[0]+=A1*S4 ; r[1]+=B1*S4 ; r[2]+=C1*S4 ; //setup.Field(r,H1);
778  pos.SetXYZ(r[0],r[1],r[2]);
779  TVector3 H1vect = GFFieldManager::getFieldVal(pos);
780  H1[0]=H1vect.X()*PS2; H1[1]=H1vect.Y()*PS2;H1[2]=H1vect.Z()*PS2; // H1 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * [(A0, B0, C0) + 2*(ax, ay, az)]
781  double A3,B3,C3,A4,B4,C4,A5,B5,C5;
782  A3 = B2*H1[2]-C2*H1[1]+A[0]; B3=C2*H1[0]-A2*H1[2]+A[1]; C3=A2*H1[1]-B2*H1[0]+A[2]; // (A2, B2, C2) x H1 + (ax, ay, az)
783  A4 = B3*H1[2]-C3*H1[1]+A[0]; B4=C3*H1[0]-A3*H1[2]+A[1]; C4=A3*H1[1]-B3*H1[0]+A[2]; // (A3, B3, C3) x H1 + (ax, ay, az)
784  A5 = A4-A[0]+A4 ; B5=B4-A[1]+B4 ; C5=C4-A[2]+C4 ; // 2*(A4, B4, C4) - (ax, ay, az)
785 
786  //
787  // Last point
788  //
789  r[0]=R[0]+S*A4 ; r[1]=R[1]+S*B4 ; r[2]=R[2]+S*C4 ; //setup.Field(r,H2);
790  pos.SetXYZ(r[0],r[1],r[2]);
791  TVector3 H2vect = GFFieldManager::getFieldVal(pos);
792  H2[0]=H2vect.X()*PS2; H2[1]=H2vect.Y()*PS2;H2[2]=H2vect.Z()*PS2; // H2 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * (A4, B4, C4)
793  double A6=B5*H2[2]-C5*H2[1], B6=C5*H2[0]-A5*H2[2], C6=A5*H2[1]-B5*H2[0]; // (A5, B5, C5) x H2
794 
795  //
796  // Test approximation quality on given step and possible step reduction
797  //
798  double EST = fabs((A1+A6)-(A3+A4))+fabs((B1+B6)-(B3+B4))+fabs((C1+C6)-(C3+C4)); // EST = ||(ABC1+ABC6)-(ABC3+ABC4)||_1 = ||(axzy x H0 + ABC5 x H2) - (ABC2 x H1 + ABC3 x H1)||_1
799  if(EST>DLT) {
800  S*=0.5;
801  stopBecauseOfMaterial = false;
802  continue;
803  }
804 
805  //
806  // Derivatives of track parameters in last point
807  //
808  if(calcCov){
809  for(int i=7; i!=ND; i+=7) { // i = 7, 14, 21, 28, 35, 42, 49; ND = 56; ND1 = 49; rows of Jacobian
810 
811  double* dR = &P[i]; // dR = (dX/dpN, dY/dpN, dZ/dpN)
812  double* dA = &P[i+3]; // dA = (dAx/dpN, dAy/dpN, dAz/dpN); N = X,Y,Z,Ax,Ay,Az,q/p
813 
814  //first point
815  double dA0 = H0[ 2]*dA[1]-H0[ 1]*dA[2]; // dA0/dp }
816  double dB0 = H0[ 0]*dA[2]-H0[ 2]*dA[0]; // dB0/dp } = dA x H0
817  double dC0 = H0[ 1]*dA[0]-H0[ 0]*dA[1]; // dC0/dp }
818 
819  if(i==ND1) {dA0+=A0; dB0+=B0; dC0+=C0;} // if last row: (dA0, dB0, dC0) := (dA0, dB0, dC0) + (A0, B0, C0)
820 
821  double dA2 = dA0+dA[0]; // }
822  double dB2 = dB0+dA[1]; // } = (dA0, dB0, dC0) + dA
823  double dC2 = dC0+dA[2]; // }
824 
825  //second point
826  double dA3 = dA[0]+dB2*H1[2]-dC2*H1[1]; // dA3/dp }
827  double dB3 = dA[1]+dC2*H1[0]-dA2*H1[2]; // dB3/dp } = dA + (dA2, dB2, dC2) x H1
828  double dC3 = dA[2]+dA2*H1[1]-dB2*H1[0]; // dC3/dp }
829 
830  if(i==ND1) {dA3+=A3-A[0]; dB3+=B3-A[1]; dC3+=C3-A[2];} // if last row: (dA3, dB3, dC3) := (dA3, dB3, dC3) + (A3, B3, C3) - (ax, ay, az)
831 
832  double dA4 = dA[0]+dB3*H1[2]-dC3*H1[1]; // dA4/dp }
833  double dB4 = dA[1]+dC3*H1[0]-dA3*H1[2]; // dB4/dp } = dA + (dA3, dB3, dC3) x H1
834  double dC4 = dA[2]+dA3*H1[1]-dB3*H1[0]; // dC4/dp }
835 
836  if(i==ND1) {dA4+=A4-A[0]; dB4+=B4-A[1]; dC4+=C4-A[2];} // if last row: (dA4, dB4, dC4) := (dA4, dB4, dC4) + (A4, B4, C4) - (ax, ay, az)
837 
838  //last point
839  double dA5 = dA4+dA4-dA[0]; // }
840  double dB5 = dB4+dB4-dA[1]; // } = 2*(dA4, dB4, dC4) - dA
841  double dC5 = dC4+dC4-dA[2]; // }
842 
843  double dA6 = dB5*H2[2]-dC5*H2[1]; // dA6/dp }
844  double dB6 = dC5*H2[0]-dA5*H2[2]; // dB6/dp } = (dA5, dB5, dC5) x H2
845  double dC6 = dA5*H2[1]-dB5*H2[0]; // dC6/dp }
846 
847  if(i==ND1) {dA6+=A6; dB6+=B6; dC6+=C6;} // if last row: (dA6, dB6, dC6) := (dA6, dB6, dC6) + (A6, B6, C6)
848 
849  dR[0]+=(dA2+dA3+dA4)*S3; dA[0] = (dA0+dA3+dA3+dA5+dA6)*P3; // dR := dR + S3*[(dA2, dB2, dC2) + (dA3, dB3, dC3) + (dA4, dB4, dC4)]
850  dR[1]+=(dB2+dB3+dB4)*S3; dA[1] = (dB0+dB3+dB3+dB5+dB6)*P3; // dA := 1/3*[(dA0, dB0, dC0) + 2*(dA3, dB3, dC3) + (dA5, dB5, dC5) + (dA6, dB6, dC6)]
851  dR[2]+=(dC2+dC3+dC4)*S3; dA[2] = (dC0+dC3+dC3+dC5+dC6)*P3;
852  }
853  }
854 
855  Way2 += S; // add stepsize to way (signed)
856  if((Way+=fabs(S))>Wmax){
857  std::cerr<<"RKTrackRep::RKutta ==> Trajectory is longer than length limit : "<<Way<<" cm !"
858  << " p/q = "<<1./P[6]<< " GeV"<<std::endl;
859  return(false);
860  }
861 
862  //
863  // Track parameters in last point
864  //
865  R[0]+=(A2+A3+A4)*S3; A[0]+=(SA[0]=(A0+A3+A3+A5+A6)*P3-A[0]); // R = R0 + S3*[(A2, B2, C2) + (A3, B3, C3) + (A4, B4, C4)]
866  R[1]+=(B2+B3+B4)*S3; A[1]+=(SA[1]=(B0+B3+B3+B5+B6)*P3-A[1]); // A = 1/3*[(A0, B0, C0) + 2*(A3, B3, C3) + (A5, B5, C5) + (A6, B6, C6)]
867  R[2]+=(C2+C3+C4)*S3; A[2]+=(SA[2]=(C0+C3+C3+C5+C6)*P3-A[2]); // SA = A_new - A_old
868  Sl=S; // last S used
869 
870  // if extrapolation has changed direction, delete the last point, because it is
871  // not a consecutive point to be used for material estimations
872  if(Ssign*S<0.) {
873  pointPaths.at(pointPaths.size()-1)+=S;
874  points.erase(points.end());
875  }
876  else{
877  pointPaths.push_back(S);
878  }
879 
880  points.push_back(TVector3(R[0],R[1],R[2]));
881 
882  double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]); // 1/|A|
883  A[0]*=CBA; A[1]*=CBA; A[2]*=CBA; // normalize A
884 
885  // Step estimation until surface and test conditions for stop of propogation
886  if(fabs(Way2)>Wmax) {
887  Dis=0.;
888  Dist=0.;
889  S=0;
890  Step=0.;
891  break;
892  }
893 
894 
895  An=A[0]*SU[0]+A[1]*SU[1]+A[2]*SU[2];
896 
897  if(fabs(An) < 1.E-6) {
898  error=true;
899  Step=0;
900  break;
901  }
902 
903  if( plane.inActive(TVector3(R[0],R[1],R[2]),TVector3(A[0],A[1],A[2]))) {
904  Dis=SU[3]-R[0]*SU[0]-R[1]*SU[1]-R[2]*SU[2];
905  Step=Dis/An;
906  }
907  else{
908  if( (O.X()-R[0])*A[0] + (O.Y()-R[1])*A[1] + (O.Z()-R[2])*A[2] >0 ){
909  Dis = sqrt((R[0]-O.X())*(R[0]-O.X())+
910  (R[1]-O.Y())*(R[1]-O.Y())+
911  (R[2]-O.Z())*(R[2]-O.Z()));
912  }
913  else{
914  Dis = -1.*sqrt((R[0]-O.X())*(R[0]-O.X())+
915  (R[1]-O.Y())*(R[1]-O.Y())+
916  (R[2]-O.Z())*(R[2]-O.Z()));
917  }
918  Step = Dis; // signed distance to surface
919  }
920 
921  if (Dis*Dist>0 && fabs(Dis)>fabs(Dist)) { // did not get closer to surface
922  error=true;
923  Step=0;
924  break;
925  }
926  Dist=Dis;
927 
928  //
929  // reset & check step size
930  //
931  // reset S to Step if extrapolation too long or in wrong direction
932  if (S*Step<0. || fabs(S)>fabs(Step)) S=Step;
933  else if (EST<DLT32 && fabs(2.*S)<=Smax) S*=2.;
934 
935  } //end of main loop
936 
937  //
938  // Output information preparation for main track parameteres
939  //
940 
941  if (!stopBecauseOfMaterial) { // linear extrapolation to surface
942  if(fabs(Sl) > 1.E-12) Sl=1./Sl; // Sl = inverted last Stepsize Sl
943  A [0]+=(SA[0]*=Sl)*Step; // Step = distance to surface
944  A [1]+=(SA[1]*=Sl)*Step; // SA*Sl = delta A / delta way; local derivative of A with respect to the length of the way
945  A [2]+=(SA[2]*=Sl)*Step; // A = A + Step * SA*Sl
946 
947  P[0] = R[0]+Step*(A[0]-.5*Step*SA[0]); // P = R + Step*(A - 1/2*Step*SA); approximation for final point on surface
948  P[1] = R[1]+Step*(A[1]-.5*Step*SA[1]);
949  P[2] = R[2]+Step*(A[2]-.5*Step*SA[2]);
950 
951  points.push_back(TVector3(P[0],P[1],P[2]));
952  pointPaths.push_back(Step);
953  }
954 
955  double CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]);
956 
957  P[3] = A[0]*CBA; // normalize A
958  P[4] = A[1]*CBA;
959  P[5] = A[2]*CBA;
960 
961  //
962  // Output derivatives of track parameters preparation
963  //
964  An = A[0]*SU[0]+A[1]*SU[1]+A[2]*SU[2];
965  fabs(An) < 1.E-6 ? An=1./An : An = 0; // 1/A_normal
966 
967  if(calcCov && !stopBecauseOfMaterial){
968  for(int i=7; i!=ND; i+=7) {
969  double* dR = &P[i]; double* dA = &P[i+3];
970  S = (dR[0]*SU[0]+dR[1]*SU[1]+dR[2]*SU[2])*An; // dR_normal / A_normal
971  dR[0]-=S*A [0]; dR[1]-=S*A [1]; dR[2]-=S*A [2];
972  dA[0]-=S*SA[0]; dA[1]-=S*SA[1]; dA[2]-=S*SA[2];
973  }
974  }
975 
976  if(error){
977  std::cerr << "RKTrackRep::RKutta ==> Do not get closer. Path = " << Way << " cm" << " p/q = " << 1./P[6] << " GeV" << std::endl;
978  return(false);
979  }
980 
981  // calculate covered distance
982  if (!stopBecauseOfMaterial) coveredDistance=Way2+Step;
983  else coveredDistance=Way2;
984 
985  return(true);
986 }
987 
988 
989 
990 
991 double RKTrackRep::Extrap( const GFDetPlane& plane, TMatrixT<double>* state, TMatrixT<double>* cov) const {
992 
993  static const int maxNumIt(2000);
994  int numIt(0);
995  bool calcCov(true);
996  if(cov==NULL) calcCov=false;
997 
998  double *P;
999  if(calcCov) {P = new double[56]; memset(P,0x00,56*sizeof(double));}
1000  else {P = new double[7];} // not needed memset(P,0x00,7*sizeof(double));};
1001 
1002  for(int i=0;i<7;++i){
1003  P[i] = (*state)[i][0];
1004  }
1005 
1006  TMatrixT<double> jac(7,7);
1007  TMatrixT<double> jacT(7,7);
1008  TMatrixT<double> oldCov(7,7);
1009  if(calcCov) oldCov=(*cov);
1010  double coveredDistance(0.);
1011  double sumDistance(0.);
1012 
1013  while(true){
1014  if(numIt++ > maxNumIt){
1015  GFException exc("RKTrackRep::Extrap ==> maximum number of iterations exceeded",
1016  __LINE__,__FILE__);
1017  exc.setFatal();
1018  delete[] P;
1019  throw exc;
1020  }
1021 
1022  if(calcCov){
1023  memset(&P[7],0x00,49*sizeof(double));
1024  for(int i=0; i<6; ++i){
1025  P[(i+1)*7+i] = 1.;
1026  }
1027  P[55] = (*state)[6][0];
1028  }
1029 
1030  //double dir(1.); //[R.K. 01/2017] unused variable?
1031  {
1032  TVector3 Pvect(P[0],P[1],P[2]); //position
1033  TVector3 Avect(P[3],P[4],P[5]); //direction
1034  TVector3 dist = plane.dist(Pvect); //from point to plane
1035  //if(dist*Avect<0.) dir=-1.; //[R.K. 01/2017] unused variable?
1036  }
1037 
1038  TVector3 directionBefore(P[3],P[4],P[5]); // direction before propagation
1039  directionBefore.SetMag(1.);
1040 
1041  // propagation
1042  std::vector<TVector3> points;
1043  std::vector<double> pointPaths;
1044  if( ! this->RKutta(plane, P, coveredDistance, points, pointPaths, -1., calcCov) ) { // maxLen currently not used
1045  GFException exc("RKTrackRep::Extrap ==> Runge Kutta propagation failed",__LINE__,__FILE__);
1046  exc.setFatal(); // stops propagation; faster, but some hits will be lost
1047  delete[] P;
1048  throw exc;
1049  }
1050 
1051  TVector3 directionAfter(P[3],P[4],P[5]); // direction after propagation
1052  directionAfter.SetMag(1.);
1053 
1054  sumDistance+=coveredDistance;
1055 
1056  // filter Points
1057  std::vector<TVector3> pointsFilt(1, points.at(0));
1058  std::vector<double> pointPathsFilt(1, 0.);
1059  // only if in right direction
1060  for(unsigned int i=1;i<points.size();++i){
1061  if (pointPaths.at(i) * coveredDistance > 0.) {
1062  pointsFilt.push_back(points.at(i));
1063  pointPathsFilt.push_back(pointPaths.at(i));
1064  }
1065  else {
1066  pointsFilt.back() = points.at(i);
1067  pointPathsFilt.back() += pointPaths.at(i);
1068  }
1069  // clean up tiny steps
1070  int position = pointsFilt.size()-1; // position starts with 0
1071  if (fabs(pointPathsFilt.back()) < MINSTEP && position > 1) {
1072  pointsFilt.at(position-1) = pointsFilt.at(position);
1073  pointsFilt.pop_back();
1074  pointPathsFilt.at(position-1) += pointPathsFilt.at(position);
1075  pointPathsFilt.pop_back();
1076  }
1077  }
1078 
1079  //consistency check
1080  double checkSum(0.);
1081  for(unsigned int i=0;i<pointPathsFilt.size();++i){
1082  checkSum+=pointPathsFilt.at(i);
1083  }
1084  if(fabs(checkSum-coveredDistance)>1.E-7){
1085  GFException exc("RKTrackRep::Extrap ==> fabs(checkSum-coveredDistance)>1.E-7",__LINE__,__FILE__);
1086  exc.setFatal();
1087  delete[] P;
1088  throw exc;
1089  }
1090 
1091  if(calcCov){ //calculate Jacobian jac
1092  for(int i=0;i<7;++i){
1093  for(int j=0;j<7;++j){
1094  if(i<6) jac[i][j] = P[ (i+1)*7+j ];
1095  else jac[i][j] = P[ (i+1)*7+j ]/P[6];
1096  }
1097  }
1098  jacT = jac;
1099  jacT.T();
1100  }
1101 
1102  TMatrixT<double> noise(7,7); // zero everywhere by default
1103 
1104  // call MatEffects
1105  double momLoss; // momLoss has a sign - negative loss means momentum gain
1106 
1107  momLoss = GFMaterialEffects::getInstance()->effects(pointsFilt,
1108  pointPathsFilt,
1109  fabs(fCharge/P[6]), // momentum
1110  fPdg,
1111  calcCov,
1112  &noise,
1113  &jac,
1114  &directionBefore,
1115  &directionAfter);
1116 
1117  if(fabs(P[6])>1.E-10){ // do momLoss only for defined 1/momentum .ne.0
1118  P[6] = fCharge/(fabs(fCharge/P[6])-momLoss);
1119  }
1120 
1121  if(calcCov){ //propagate cov and add noise
1122  if(!(oldCov < 1.E200)){
1123  GFException exc("RKTrackRep::Extrap ==> covariance matrix exceeds numerical limits",__LINE__,__FILE__);
1124  exc.setFatal();
1125  delete[] P;
1126  throw exc;
1127  }
1128  oldCov = *cov;
1129  *cov = jacT*((oldCov)*jac)+noise;
1130  }
1131 
1132 
1133  //we arrived at the destination plane, if we point to the active area
1134  //of the plane (if it is finite), and the distance is below threshold
1135  if( plane.inActive(TVector3(P[0],P[1],P[2]),TVector3(P[3],P[4],P[5]))) {
1136  if(plane.distance(P[0],P[1],P[2])<MINSTEP) break;
1137  }
1138  }
1139  (*state)[0][0] = P[0]; (*state)[1][0] = P[1];
1140  (*state)[2][0] = P[2]; (*state)[3][0] = P[3];
1141  (*state)[4][0] = P[4]; (*state)[5][0] = P[5];
1142  (*state)[6][0] = P[6];
1143 
1144  delete[] P;
1145  return sumDistance;
1146 }
1147 
1148 void RKTrackRep::getPosMomCov(const GFDetPlane& pl, TVector3& pos, TVector3& mom, TMatrixT<double>& cov){
1149 
1150  const TVector3 U=pl.getU();
1151  const TVector3 V=pl.getV();
1152  const TVector3 W=U.Cross(V);
1153 
1154  TMatrixT<double> statePred=fState;
1155  TMatrixT<double> covPred=fCov;
1156  double spu=fSpu;
1157  if(pl!=fRefPlane)
1158  {
1159  extrapolate(pl, statePred, covPred);
1160  spu=fCacheSpu;
1161  }
1162 
1163  const double p=fCharge/statePred[0][0]; //Magnitude of the momentum.
1164  const TVector3 pTilde=spu*(W+statePred[1][0]*U+statePred[2][0]*V);
1165  const double pTildeMagnitude=pTilde.Mag();
1166 
1167  TMatrixT<double> jacobian(6, 5); //Jacobian matrix \frac{\partial\left(x, y, z, p_x, p_y, p_z\right)}{\partial\left(\frac{q}{p}, u^\prime, v^\prime, u, v\right)}
1168  jacobian[0][3]=U.X();
1169  jacobian[0][4]=V.X();
1170  jacobian[1][3]=U.Y();
1171  jacobian[1][4]=V.Y();
1172  jacobian[2][3]=U.Z();
1173  jacobian[2][4]=V.Z();
1174  jacobian[3][1]=p*spu/pTildeMagnitude*(U.X()-pTilde.X()*(U*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1175  jacobian[3][2]=p*spu/pTildeMagnitude*(V.X()-pTilde.X()*(V*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1176  jacobian[4][1]=p*spu/pTildeMagnitude*(U.Y()-pTilde.Y()*(U*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1177  jacobian[4][2]=p*spu/pTildeMagnitude*(V.Y()-pTilde.Y()*(V*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1178  jacobian[5][1]=p*spu/pTildeMagnitude*(U.Z()-pTilde.Z()*(U*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1179  jacobian[5][2]=p*spu/pTildeMagnitude*(V.Z()-pTilde.Z()*(V*pTilde)/(pTildeMagnitude*pTildeMagnitude));
1180 
1181  TMatrixT<double> transposedJacobian=jacobian;
1182  transposedJacobian.T();
1183 
1184  cov.ResizeTo(6, 6);
1185 
1186  pos=pl.getO()+statePred[3][0]*U+statePred[4][0]*V;
1187  mom=p/pTildeMagnitude*pTilde;
1188  cov=jacobian*covPred*transposedJacobian;
1189 }
1190 
1191 
1192 double RKTrackRep::stepalong(double h, TVector3& pos, TVector3& dir){
1193 
1194  TVector3 dest;
1195 
1196  static const int maxIt(30);
1197 
1198  TVector3 o=fRefPlane.getO();
1199  TVector3 u=fRefPlane.getU();
1200  TVector3 v=fRefPlane.getV();
1201  TVector3 w=u.Cross(v);
1202 
1203  TVector3 pTilde = fSpu* (w + fState[1][0] * u + fState[2][0] * v);
1204  pTilde.SetMag(1.);
1205 
1206  TVector3 point = o + fState[3][0]*u + fState[4][0]*v;
1207 
1208  TMatrixT<double> state7(7,1);
1209  state7[0][0] = point.X();
1210  state7[1][0] = point.Y();
1211  state7[2][0] = point.Z();
1212  state7[3][0] = pTilde.X();
1213  state7[4][0] = pTilde.Y();
1214  state7[5][0] = pTilde.Z();
1215  state7[6][0] = fState[0][0];
1216 
1217  double coveredDistance(0.);
1218 
1219  GFDetPlane pl;
1220  int iterations(0);
1221 
1222  while(true){
1223  pos.SetXYZ(state7[0][0], state7[1][0], state7[2][0]);
1224  dir.SetXYZ(state7[3][0], state7[4][0], state7[5][0]);
1225  dir.SetMag(1.);
1226 
1227  dest = pos + (h - coveredDistance) * dir;
1228 
1229  pl.setON(dest, dir);
1230  coveredDistance += this->Extrap(pl,&state7);
1231 
1232  if(fabs(h - coveredDistance)<MINSTEP) break;
1233  if(++iterations == maxIt) {
1234  GFException exc("RKTrackRep::stepalong ==> maximum number of iterations reached",__LINE__,__FILE__);
1235  throw exc;
1236  }
1237  }
1238 
1239  pos.SetXYZ(state7[0][0], state7[1][0], state7[2][0]);
1240  dir.SetXYZ(state7[3][0], state7[4][0], state7[5][0]);
1241  dir.SetMag(1.);
1242 
1243  return coveredDistance;
1244 
1245 }
1246 
1247 
1248 
1249 ClassImp(RKTrackRep)
GFAbsTrackRep::getMom
TVector3 getMom()
Definition: GFAbsTrackRep.h:221
pos
TVector3 pos
Definition: GammaSpectraAnalysis.C:43
GFAbsTrackRep
Base Class for genfit track representations. Defines interface for track parameterizations.
Definition: GFAbsTrackRep.h:80
RKTrackRep::setData
void setData(const TMatrixT< double > &st, const GFDetPlane &pl, const TMatrixT< double > *cov=NULL, const TMatrixT< double > *aux=NULL)
Sets state, plane and (optionally) covariance.
Definition: RKTrackRep.cxx:43
GFDetPlane::getV
TVector3 getV() const
Definition: GFDetPlane.h:77
r
double r
Definition: RiemannTest.C:14
RKTrackRep::Extrap
double Extrap(const GFDetPlane &plane, TMatrixT< double > *state, TMatrixT< double > *cov=NULL) const
Handles propagation and material effects.
Definition: RKTrackRep.cxx:991
i
Int_t i
Definition: run_full.C:25
GFDetPlane
Detector plane genfit geometry class.
Definition: GFDetPlane.h:59
RKTrackRep::getPosMom
void getPosMom(const GFDetPlane &, TVector3 &pos, TVector3 &mom)
Gets position and momentum in the plane.
Definition: RKTrackRep.cxx:312
GFDetPlane::setU
void setU(const TVector3 &u)
Definition: GFDetPlane.cxx:105
exit
exit(0)
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
s
TLorentzVector s
Definition: Pnd2DStar.C:50
GFDetPlane::setON
void setON(const TVector3 &o, const TVector3 &n)
Definition: GFDetPlane.cxx:145
GFDetPlane::setO
void setO(const TVector3 &o)
Definition: GFDetPlane.cxx:92
GFDetPlane::distance
double distance(TVector3 &) const
Definition: GFDetPlane.cxx:343
GFFieldManager::getFieldVal
static TVector3 getFieldVal(const TVector3 &x)
Definition: GFFieldManager.h:52
RKTrackRep::extrapolateToLine
void extrapolateToLine(const TVector3 &point1, const TVector3 &point2, TVector3 &poca, TVector3 &dirInPoca, TVector3 &poca_onwire)
This method extrapolates to the point of closest approach to a line.
Definition: RKTrackRep.cxx:390
mom
Double_t mom
Definition: plot_dirc.C:14
Y
double Y
Definition: anaLmdDigi.C:68
v
__m128 v
Definition: P4_F32vec4.h:4
error
static void error(int no)
Definition: ranlxd.cxx:419
RKTrackRep::RKutta
bool RKutta(const GFDetPlane &plane, double *P, double &coveredDistance, std::vector< TVector3 > &points, std::vector< double > &pointLengths, const double &maxLen=-1, bool calcCov=true) const
Propagates the particle through the magnetic field.
Definition: RKTrackRep.cxx:640
p
Double_t p
Definition: anasim.C:58
GFDetPlane::setNormal
void setNormal(TVector3 n)
Definition: GFDetPlane.cxx:156
GFTrackCand::getPosSeed
TVector3 getPosSeed() const
get the seed value for track: pos
Definition: GFTrackCand.h:131
RKTrackRep::~RKTrackRep
virtual ~RKTrackRep()
Definition: RKTrackRep.cxx:69
GFTrackCand::getDirSeed
TVector3 getDirSeed() const
get the seed value for track: direction
Definition: GFTrackCand.h:133
W
#define W
Definition: createSTT.C:76
RKTrackRep
Track Representation module based on a Runge-Kutta algorithm including a full material model...
Definition: tracking/GenfitTools/trackrep/RKTrackRep/RKTrackRep.h:49
RKTrackRep::getPosMomCov
void getPosMomCov(const GFDetPlane &pl, TVector3 &pos, TVector3 &mom, TMatrixT< double > &cov)
method which gets position, momentum and 6x6 covariance matrix
Definition: RKTrackRep.cxx:1148
point
TClonesArray * point
Definition: anaLmdDigi.C:29
GFTrackCand::getDirError
TVector3 getDirError() const
get the seed value for track: error on direction (standard deviation)
Definition: GFTrackCand.h:139
GFAbsTrackRep::fRefPlane
GFDetPlane fRefPlane
Definition: GFAbsTrackRep.h:112
GFTrackCand
Track candidate – a list of cluster indices.
Definition: GFTrackCand.h:55
MINSTEP
#define MINSTEP
Definition: RKTrackRep.cxx:40
GFAbsTrackRep::getPos
TVector3 getPos()
Definition: GFAbsTrackRep.h:220
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
GFTrackCand::getPosError
TVector3 getPosError() const
get the seed value for track: error on pos (standard deviation)
Definition: GFTrackCand.h:137
RKTrackRep::stepalong
double stepalong(double h, TVector3 &point, TVector3 &dir)
make step of h cm along the track, returns the tracklength spanned in this extrapolation ...
Definition: RKTrackRep.cxx:1192
RKTrackRep::extrapolateToPoint
void extrapolateToPoint(const TVector3 &pos, TVector3 &poca, TVector3 &dirInPoca)
This method is to extrapolate the track to point of closest approach to a point in space...
Definition: RKTrackRep.cxx:329
GFTrackCand::getPdgCode
int getPdgCode() const
get the PDG code
Definition: GFTrackCand.h:141
GFMaterialEffects::stepper
double stepper(const double &maxDist, const double &posx, const double &posy, const double &posz, const double &dirx, const double &diry, const double &dirz, const double &mom, const int &pdg)
Returns maximum length so that a specified momentum loss will not be exceeded.
Definition: GFMaterialEffects.cxx:139
X
double X
Definition: anaLmdDigi.C:68
GFDetPlane::setV
void setV(const TVector3 &v)
Definition: GFDetPlane.cxx:118
RKTrackRep::poca2Line
TVector3 poca2Line(const TVector3 &extr1, const TVector3 &extr2, const TVector3 &point) const
Definition: RKTrackRep.cxx:376
P
GeV c P
Definition: ReadMesDecCorrRes_PiT.C:42
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:45
GFDetPlane::getU
TVector3 getU() const
Definition: GFDetPlane.h:76
ClassImp
ClassImp(PndAnaContFact)
Z
double Z
Definition: anaLmdDigi.C:68
GFException::setFatal
void setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:77
t
TTree * t
Definition: bump_analys.C:13
GFDetPlane::inActive
bool inActive(const TVector3 &point, const TVector3 &dir) const
intersect in the active area? C.f. GFAbsFinitePlane
Definition: GFDetPlane.h:132
GFAbsTrackRep::fState
TMatrixT< double > fState
The vector of track parameters.
Definition: GFAbsTrackRep.h:88
GFAbsTrackRep::setData
virtual void setData(const TMatrixT< double > &st, const GFDetPlane &pl, const TMatrixT< double > *cov=NULL, const TMatrixT< double > *aux=NULL)
Puts the track representation in a given state.
Definition: GFAbsTrackRep.h:267
RKTrackRep::setPDG
void setPDG(int)
Set PDG particle code.
Definition: RKTrackRep.cxx:274
GFMaterialEffects::getInstance
static GFMaterialEffects * getInstance()
Definition: GFMaterialEffects.cxx:63
GFDetPlane::getNormal
TVector3 getNormal() const
Definition: GFDetPlane.cxx:138
RKTrackRep::extrapolate
double extrapolate(const GFDetPlane &, TMatrixT< double > &statePred, TMatrixT< double > &covPred)
returns the tracklength spanned in this extrapolation
Definition: RKTrackRep.cxx:442
noise
double noise
Definition: outdated/dc1/stt/runreco.C:105
R
Double_t R
Definition: checkhelixhit.C:61
GFMaterialEffects::effects
double effects(const std::vector< TVector3 > &points, const std::vector< double > &pointPaths, const double &mom, const int &pdg, const bool &doNoise=false, TMatrixT< double > *noise=NULL, const TMatrixT< double > *jacobian=NULL, const TVector3 *directionBefore=NULL, const TVector3 *directionAfter=NULL)
Calculates energy loss in the travelled path, optional calculation of noise matrix.
Definition: GFMaterialEffects.cxx:76
GFDetPlane::getO
TVector3 getO() const
Definition: GFDetPlane.h:75
RKTrackRep::getAuxInfo
const TMatrixT< double > * getAuxInfo(const GFDetPlane &pl)
Get auxillary information from the track representation.
Definition: RKTrackRep.cxx:57
GFAbsTrackRep::fCov
TMatrixT< double > fCov
The covariance matrix.
Definition: GFAbsTrackRep.h:91
GFDetPlane::dist
TVector3 dist(const TVector3 &point) const
Definition: GFDetPlane.cxx:207
Generated on Sun Aug 18 2019 03:01:58 for FairRoot/PandaRoot by   doxygen 1.8.8