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