EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GeaneTrackRep.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GeaneTrackRep.cxx
1 //-----------------------------------------------------------
2 // File and Version Information:
3 // $Id$
4 //
5 // Description:
6 // Implementation of class GeaneTrackRep
7 // see GeaneTrackRep.hh for details
8 //
9 // Environment:
10 // Software developed for the PANDA Detector at FAIR.
11 //
12 // Author List:
13 // Sebastian Neubert TUM (original author)
14 //
15 //
16 //-----------------------------------------------------------
17 
18 // Panda Headers ----------------------
19 
20 // This Class' Header ------------------
21 #include "GeaneTrackRep.h"
22 #include "FairGeaneUtil.h"
23 #include "FairTrackParH.h"
24 
25 // C/C++ Headers ----------------------
26 #include <iostream>
27 #include <cmath>
28 
29 // Collaborating Class Headers --------
30 #include "GFAbsRecoHit.h"
31 #include "GFException.h"
32 #include "FairGeanePro.h"
33 
34 // Class Member definitions -----------
35 
36 #define THETACUT 0.4
37 #define EPSILON 1E-4
38 
39 GeaneTrackRep::GeaneTrackRep()
40  : GFAbsTrackRep(5), _pdg(211),_backw(0), _spu(1)
41 {
42 
43 }
44 
45 GeaneTrackRep::GeaneTrackRep(FairGeanePro* geane,
46  const GFDetPlane& plane,
47  const TVector3& mom,
48  const TVector3& poserr,
49  const TVector3& momerr,
50  double q,
51  int PDGCode)
52  : GFAbsTrackRep(5), _geane(geane), _pdg(PDGCode), _backw(0), _spu(1)
53 {
54  FairTrackParP par(plane.getO(),mom,poserr,momerr,(int)TMath::Sign(1.0, q),plane.getO(),plane.getU(),plane.getV());
55 
56  _spu=par.GetSPU(); // direction of the momentum
57 
58  fState[0][0]=par.GetQp();
59  fState[1][0]=par.GetTV();
60  fState[2][0]=par.GetTW();
61  fState[3][0]=par.GetV();
62  fState[4][0]=par.GetW();
63 
64  // blow up cov-array: ROOT does not support init with symmetric data
65  // See ROOT docu source-file for TMatrixTSym
66  // i=row, j=collumn
67  double* covarray=par.GetCov();
68  int count=0;
69  for(int i=0;i<5;++i){
70  for(int j=i;j<5;++j){
71  fCov[i][j]=covarray[count];
72  if(i!=j)fCov[j][i]=covarray[count];
73  ++count;
74  }
75  }
76  fRefPlane=plane;
77 }
78 
79 GeaneTrackRep::GeaneTrackRep(FairGeanePro* geane,
80  const GFDetPlane& plane,
81  const TVector3& mom,
82  const TVector3& poserr,
83  const TVector3& momerr,
84  int q,
85  int PDGCode)
86  : GFAbsTrackRep(5), _geane(geane), _pdg(PDGCode), _backw(0), _spu(1)
87 {
88  FairTrackParP par(plane.getO(),mom,poserr,momerr,q,plane.getO(),plane.getU(),plane.getV());
89 
90  _spu=par.GetSPU(); // direction of the momentum
91 
92  fState[0][0]=par.GetQp();
93  fState[1][0]=par.GetTV();
94  fState[2][0]=par.GetTW();
95  fState[3][0]=par.GetV();
96  fState[4][0]=par.GetW();
97 
98  // blow up cov-array: ROOT does not support init with symmetric data
99  // See ROOT docu source-file for TMatrixTSym
100  // i=row, j=collumn
101  double* covarray=par.GetCov();
102  int count=0;
103  for(int i=0;i<5;++i){
104  for(int j=i;j<5;++j){
105  fCov[i][j]=covarray[count];
106  if(i!=j)fCov[j][i]=covarray[count];
107  ++count;
108  }
109  }
110  fRefPlane=plane;
111 }
112 
113 // GeaneTrackRep::GeaneTrackRep(const GeaneTrackRep& rep)
114 // : GFAbsTrackRep(rep)
115 // {
116 // _geane=rep._geane;
117 // }
118 
119 
120 GeaneTrackRep::~GeaneTrackRep()
121 {
122 
123 }
124 
125 
126 
127 
128 double
129 GeaneTrackRep::extrapolate(const GFDetPlane& pl,
130  TMatrixT<double>& statePred)
131 {
132  TMatrixT<double> covPred(5,5);
133  return extrapolate(pl,statePred,covPred);
135 }
136 
137 
138 double
139 GeaneTrackRep::extrapolate(const GFDetPlane& pl,
140  TMatrixT<double>& statePred,
141  TMatrixT<double>& covPred)
142 {
143  //printf("@@@ Geane-1\n");
144  if(fabs(getMom(fRefPlane).Theta()/TMath::Pi()*180.) < THETACUT){
145  GFException exc("GEANE propagation not possible for p.theta<THETACUT",__LINE__,__FILE__);
146  exc.setFatal();
147  throw exc;
148  }
149  statePred.ResizeTo(fDimension,1);
150  covPred.ResizeTo(fDimension,fDimension);
151  TVector3 o=pl.getO();
152  TVector3 u=pl.getU();
153  TVector3 v=pl.getV();
154 
155  TVector3 ofrom=fRefPlane.getO();
156  TVector3 ufrom=fRefPlane.getU();
157  TVector3 vfrom=fRefPlane.getV();
158 
159  _geane->PropagateFromPlane(ufrom,vfrom);
160  _geane->PropagateToPlane(o,u,v);
161 
162  FairTrackParP result;
163  FairTrackParH result2;
164 
165  //std::cout<<"Before prop:"<<std::endl;
166  //Print();
167 
168 
169 
170  double cova[15];
171  int count=0;;
172  for(int i=0; i<5;++i){
173  for(int j=i;j<5;++j){
174  cova[count++]=fCov[i][j];
175  }
176  }
177  //protect against low momentum:
178  if(fabs(fState[0][0])>10){
179  GFException exc("GeaneTrackRep: PROTECT AGAINST LOW MOMENTA",__LINE__,__FILE__);
180  exc.setFatal();
181  throw exc;
182  }
183 
184  //printf("@@@ Geane-2\n");
185  checkState();
186 
187  //printf("@@@ Geane-3\n");
188 
189  FairTrackParP par(fState[3][0],fState[4][0],fState[1][0],fState[2][0],fState[0][0],cova,ofrom,ufrom,vfrom,_spu);
190  bool backprop=_backw<0;
191  if(_backw==0){
192  //Try to guess if we are doing a forward or backward step:
193  TVector3 pos(par.GetX(),par.GetY(),par.GetZ());
194  TVector3 dir=pl.dist(pos); // direction from pos to plane;
195  //Assume B=(0,0,BZ) -> compare signs of dir.Z and mom.Z:
196  //backprop= (dir.Z()*par.GetPz())<0 ? true : false;
197  backprop= (dir*getMom(fRefPlane))<0;
198  }
199  if(backprop){
200  _geane->setBackProp();
201  //std::cout<<"GEANETRACKREP: USING BACKPROPAGATION!" << std::endl;
202  }
203 
204  //printf("@@@ Geane-3a %d %d\n", _backw, backprop);
205  Bool_t prop = kTRUE;
206  prop = _geane->Propagate(&par,&result,_pdg); //211
207  //printf("@@@ Geane-3b: %d\n", prop);
208  if (prop==kFALSE){
209  GFException exc("GEANE propagation failed",__LINE__,__FILE__);
210  //exc.setFatal();
211  throw exc;
212  }
213 
214  //printf("@@@ Geane-4\n");
215 
216  double l=_geane->GetLengthAtPCA();
217 
218  statePred[0][0]=result.GetQp();
219  statePred[1][0]=result.GetTV();
220  statePred[2][0]=result.GetTW();
221  statePred[3][0]=result.GetV();
222  statePred[4][0]=result.GetW();
223 
224 
225 
226  double* rescov=result.GetCov();
227  count=0;
228  for(int i=0;i<5;++i){
229  for(int j=i;j<5;++j){
230  covPred[i][j]=rescov[count];
231  if(i!=j)covPred[j][i]=rescov[count];
232  ++count;
233  }
234  }
235 
236  // if(result.GetSPU()!=_spu)std::cout<<"SPU HAS CHANGED! "<<_spu<<" --> "<<result.GetSPU()<<std::endl;
237  _spu=result.GetSPU();
238 
239  //std::cout<<"AFTER EXTRAPOLATE:"<<std::endl;
240  //result.Print();
241  //pl.Print();
242  //statePred.Print();
243  //covPred.Print();
244 
245 
246 
247  return l;
248 }
249 
250 
251 
252 void
253 GeaneTrackRep::extrapolateToPoint(const TVector3& pos,
254  TVector3& poca,
255  TVector3& dirInPoca){
256  if(fabs(getMom(fRefPlane).Theta()/TMath::Pi()*180.) < THETACUT){
257  GFException exc("GEANE propagation not possible for p.theta<THETACUT",__LINE__,__FILE__);
258  exc.setFatal();
259  throw exc;
260  }
261  int dim = getDim();
262  TMatrixT<double> statePred(dim,1);
263  TMatrixT<double> covPred(dim,dim);
264  //std::cout<<"GeaneTrackRep::extrapolateToPoint"<<std::endl;
265  //fRefPlane.Print();
266 
267  TVector3 ofrom=fRefPlane.getO();
268  TVector3 ufrom=fRefPlane.getU();
269  TVector3 vfrom=fRefPlane.getV();
270 
271  _geane->SetPoint(pos);
272  _geane->PropagateFromPlane(ufrom,vfrom);
273 
274  double cova[15];
275  int count=0;;
276  for(int i=0; i<5;++i){
277  for(int j=i;j<5;++j){
278  cova[count++]=fCov[i][j];
279  }
280  }
281  //protect against low momentum:
282  if(fabs(fState[0][0])>10){
283  GFException exc("GeaneTrackRep: PROTECT AGAINST LOW MOMENTA",__LINE__,__FILE__);
284  exc.setFatal();
285  throw exc;
286  }
287 
288  checkState();
289 
290  FairTrackParP par(fState[3][0],fState[4][0],fState[1][0],fState[2][0],fState[0][0],cova,ofrom,ufrom,vfrom,_spu);
291  //par.Print();
292  bool backprop=_backw<0;
293  if(_backw==0){
294  // check if new point is after or before my position
295  TVector3 dir=fRefPlane.dist(pos); // direction from pos to plane;
296  backprop= (dir*getMom(fRefPlane))>0;
297  }
298  if(!backprop){ // point lies in same direction of flight as momentum
299  //std::cout<<" Propagate in flight direction"<<std::endl;
300  _geane->PropagateToVirtualPlaneAtPCA(1);
301  }
302  else{
303  //std::cout<<" backPropagate"<<std::endl;
304  _geane->BackTrackToVirtualPlaneAtPCA(1);
305  }
306 
307  FairTrackParP result;
308  Bool_t prop = kTRUE;
309  prop = _geane->Propagate(&par,&result,_pdg); //211
310  if (prop==kFALSE) {
311  GFException exc("GEANE propagation failed",__LINE__,__FILE__);
312  //exc.setFatal();
313  throw exc;
314  //pl=fRefPlane;
315  //return pos;
316  }
317 
318  statePred[0][0]=result.GetQp();
319  statePred[1][0]=result.GetTV();
320  statePred[2][0]=result.GetTW();
321  statePred[3][0]=result.GetV();
322  statePred[4][0]=result.GetW();
323 
324  double* rescov=result.GetCov();
325  count=0;
326  for(int i=0;i<5;++i){
327  for(int j=i;j<5;++j){
328  covPred[i][j]=rescov[count];
329  if(i!=j)covPred[j][i]=rescov[count];
330  ++count;
331  }
332  }
333 
334  poca.SetXYZ(result.GetX(),result.GetY(),result.GetZ());
335  dirInPoca = result.GetJVer().Cross( result.GetKVer() );
336 }
337 
338 
339 void
340 GeaneTrackRep::extrapolateToLine(const TVector3& point1,
341  const TVector3& point2,
342  TVector3& poca,
343  TVector3& dirInPoca,
344  TVector3& poca_onwire)
345 {
346  if(fabs(getMom(fRefPlane).Theta()/TMath::Pi()*180.) < THETACUT){
347  GFException exc("GEANE propagation not possible for p.theta<THETACUT",__LINE__,__FILE__);
348  exc.setFatal();
349  throw exc;
350  }
351 
352  // call propagation to closest approach to a wire
353  Int_t pca = 2;
354 
355  // calculate a very large track length
356  TVector3 start = getPos(fRefPlane);
357  Double_t distance1, distance2;
358  distance1 = (point1 - start).Mag();
359  distance2 = (point2 - start).Mag();
360  Double_t maxdistance;
361  if(distance1 < distance2) maxdistance = distance2;
362  else maxdistance = distance1;
363  maxdistance *= 2.;
364 
365  // variables for FindPCA:
366  TVector3 point(0,0,0);
367  Double_t Rad = 0.;
368  // poca = vpf = point of closest approach on track
369  // poca_onwire = vwi = point of closest approach on wire
370  Double_t Di = 0.;
371  Float_t trklength = 0.;
372 
373  // covariance matrix
374  FairGeaneUtil util;
375  Double_t cov55[5][5];
376  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) cov55[i][j] = fCov[i][j];
377  Double_t cova[15];
378  util.FromMat25ToVec15(cov55, cova);
379 
380  TVector3 o = fRefPlane.getO();
381  TVector3 dj = fRefPlane.getU();
382  TVector3 dk = fRefPlane.getV();
383 
384  FairTrackParP par(fState[3][0],fState[4][0],fState[1][0],fState[2][0],fState[0][0],cova,o,dj,dk,_spu);
385 
386  // get propagation direction
387  Int_t direction = getPropDir();
388 
389  _geane->ActualFindPCA(pca, &par, direction);
390  Int_t findpca = _geane->FindPCA(pca, _pdg, point, point1, point2, maxdistance, Rad, poca, poca_onwire, Di, trklength);
391 
392  if(findpca != 0) {
393  GFException exc("findpca failure", __LINE__,__FILE__);
394  throw exc;
395  }
396 
397  // dir in poca not filled now
398  dirInPoca.SetXYZ(0., 0., 0.);
399 
400 }
401 
402 
403 
404 
405 
406 
407 TVector3
408 GeaneTrackRep::getPocaOnLine(const TVector3& p1, const TVector3& p2, bool back){
409 
410  //std::cout<<"GeaneTrackRep::getPocaToWire"<<std::endl;
411 
412  TVector3 ofrom=fRefPlane.getO();
413  TVector3 ufrom=fRefPlane.getU();
414  TVector3 vfrom=fRefPlane.getV();
415 
416  _geane->SetWire(p1,p2);
417  _geane->PropagateFromPlane(ufrom,vfrom);
418  double cova[15];
419  int count=0;;
420  for(int i=0; i<5;++i){
421  for(int j=i;j<5;++j){
422  cova[count++]=fCov[i][j];
423  }
424  }
425  // protect against low momentum:
426  if(fabs(fState[0][0])>10){
427  GFException exc("GeaneTrackRep: PROTECT AGAINST LOW MOMENTA",__LINE__,__FILE__);
428  exc.setFatal();
429  throw exc;
430  }
431 
432  checkState();
433 
434  FairTrackParP par(fState[3][0],fState[4][0],fState[1][0],fState[2][0],fState[0][0],cova,ofrom,ufrom,vfrom,_spu);
435 
436 
437  if(!back){ // point lies in same direction of flight as momentum
438  //std::cout<<" Propagate in flight direction"<<std::endl;
439  _geane->PropagateToVirtualPlaneAtPCA(2); // option 2 means wire!
440  }
441  else{
442  //std::cout<<" backPropagate"<<std::endl;
443  _geane->BackTrackToVirtualPlaneAtPCA(2);
444  }
445 
446  FairTrackParP result;
447  Bool_t prop = kTRUE;
448 
449  prop = _geane->Propagate(&par,&result,_pdg);
450  if (prop==kFALSE) {
451  GFException exc("GEANE propagation failed",__LINE__,__FILE__);
452  //exc.setFatal();
453  throw exc;
454  }
455 
456  return _geane->GetPCAOnWire();
457 }
458 
459 
460 
461 
462 
463 
464 TVector3
465 GeaneTrackRep::getPos(const GFDetPlane& pl)
466 {
467  TMatrixT<double> statePred(fState);
468  if(pl!=fRefPlane)extrapolate(pl,statePred);
469  return pl.getO()+(statePred[3][0]*pl.getU())+(statePred[4][0]*pl.getV());
470 }
471 
472 TVector3
473 GeaneTrackRep::getMom(const GFDetPlane& pl)
474 {
475  TMatrixT<double> statePred(fState);
476  if(pl!=fRefPlane)extrapolate(pl,statePred);
477  double fSPU = _spu;
478  TVector3 mom = fSPU*pl.getNormal()+fSPU*statePred[1][0]*pl.getU()+fSPU*statePred[2][0]*pl.getV();
479  mom.SetMag(1./fabs(statePred[0][0]));
480  return mom;
481 }
482 void
483 GeaneTrackRep::getPosMom(const GFDetPlane& pl,TVector3& pos, TVector3& mom)
484 {
485  TMatrixT<double> statePred(fState);
486  if(pl!=fRefPlane)extrapolate(pl,statePred);
487  mom = pl.getNormal()+statePred[1][0]*pl.getU()+statePred[2][0]*pl.getV();
488 
489  mom.SetMag(1./fabs(statePred[0][0]));
490  pos = pl.getO()+(statePred[3][0]*pl.getU())+(statePred[4][0]*pl.getV());
491 }
492 
493 
494 void
495 GeaneTrackRep::getPosMomCov(const GFDetPlane& pl,TVector3& pos,TVector3& mom,TMatrixT<double>& cov){
496  cov.ResizeTo(6,6);
497 
498  TMatrixT<double> statePred(fState), covPred(fCov);
499  if(pl!=fRefPlane)extrapolate(pl, statePred, covPred);
500 
501  // position
502  pos = pl.getO()+(statePred[3][0]*pl.getU())+(statePred[4][0]*pl.getV());
503 
504  // momentum
505  double fSPU = _spu;
506  mom = fSPU*pl.getNormal()+fSPU*statePred[1][0]*pl.getU()+fSPU*statePred[2][0]*pl.getV();
507  mom.SetMag(1./fabs(statePred[0][0]));
508 
509  // covariance matrix
510  FairGeaneUtil util;
511  // covPred 5 X 5 ==> cov55[5][5]
512  double cov55[5][5];
513  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) cov55[i][j] = covPred[i][j];
514  // cov55[5][5] ==> cov15[15]
515  double cov15[15];
516  util.FromMat25ToVec15(cov55, cov15);
517 
518  FairTrackParP parPred(statePred[3][0],
519  statePred[4][0], statePred[1][0],
520  statePred[2][0], statePred[0][0],
521  cov15,
522  pl.getO(), pl.getU(), pl.getV(),
523  _spu);
524  double cov66[6][6];
525  parPred.GetMARSCov(cov66);
526  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) cov[i][j] = cov66[i][j];
527 
528 }
529 
530 
531 void
532 GeaneTrackRep::checkState(){
533  if(fabs(fState[3][0])<1.E-4)fState[3][0]=1.E-4;
534  if(fabs(fState[4][0])<1.E-4)fState[4][0]=1.E-4;
535 
536  //if (!(fState.Abs()>1.E-15) || !(fState.Abs()<1.E50)){
537  // GFException exc("fState out of numerical bounds",__LINE__,__FILE__);
538  // exc.setFatal();
539  // throw exc;
540  //}
541 }
542 
543 
544 ClassImp(GeaneTrackRep)
545