EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GFTrackProximity.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GFTrackProximity.cxx
1 /* Copyright 2008-2009, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 #include "TVector3.h"
20 #include <iostream>
21 #include "GFTrackProximity.h"
22 #include "GFTrack.h"
23 #include "GFAbsTrackRep.h"
24 using namespace std;
25 
26 double trkDist(GFAbsTrackRep* rep1, GFAbsTrackRep* rep2){
27 
28  TVector3 d=rep1->getPos()-rep2->getPos();
29  //cout<<"trkDist"<<d.Mag()<<endl;
30  return d.Mag();
31 }
32 
33 TVector3 trackProximity(GFTrack *trk1, GFTrack *trk2){
34 
35  GFAbsTrackRep* rep1=trk1->getCardinalRep();
36  GFAbsTrackRep* rep2=trk2->getCardinalRep();
37  return trackProximity(rep1,rep2);
38 }
39 
40 TVector3 trackProximity(GFAbsTrackRep* rep1, GFAbsTrackRep* rep2){
41  // TODO: make accuracy configurable!
42  // a simple newtonian search.
43  double h=0.01;
44  double m1=999999;
45  double s1=0;
46  int steps=0;
47  double oldd=999999;
48  double d=99999;
49  // TVector3 pos,mom,poserr,momerr;
50 
51 
52  //cout<<"GFTrackProximity start ########################################################################################################################"<<endl;
53 
54 
55 
56  while(oldd>d && steps<100){
57 
58  //cout<<"GFTrackProximity::steps:"<<steps<<endl;
59  TVector3 point1,dir1;
60  TVector3 point2, norm2;
61  GFDetPlane plane1,plane2;
62  rep1->stepalong(s1,point1,dir1);
63  plane1.setO(point1);
64  //cout<<"track1 ";point1.Print();
65  plane1.setNormal(dir1);
66  rep1->GFAbsTrackRep::extrapolate(plane1);
67  rep2->extrapolateToPoint(point1,point2,norm2);
68  plane2.setO(point1);
69  plane2.setNormal(norm2);
70  rep2->GFAbsTrackRep::extrapolate(plane2);
71 
72  oldd=d;
73  ++steps;
74  d=trkDist(rep1,rep2);
75  //cout<<"dist"<<d<<endl;
76  //----------------------------------------------------------------------
77  rep1->stepalong(h,point1,dir1);
78  //cout<<"track1 ";point1.Print();
79  plane1.setO(point1);
80  plane1.setNormal(dir1);
81  rep1->GFAbsTrackRep::extrapolate(plane1);
82  rep2->extrapolateToPoint(point1,point2,norm2);
83  plane2.setO(point1);
84  plane2.setNormal(norm2);
85  rep2->GFAbsTrackRep::extrapolate(plane2);
86 
87  double d1=trkDist(rep1,rep2);
88  //cout<<"dist"<<d1<<endl;
89  m1=(d1-d)/h;
90  if(TMath::Abs(m1)<1E-4)break;
91  s1=-d/m1;
92  // limit stepsize to max 180cm
93  if(TMath::Abs(s1)>200)s1= s1>0 ? 180 : -180;
94  //std:://cout<<"d="<<d<<" d1="<<d1<<" m1="<<m1<<" s1="<<s1<<std::endl;
95 
96  }
97  TVector3 point1,dir1;
98  GFDetPlane plane1, plane2;
99  TVector3 point2, norm2;
100  rep1->stepalong(s1,point1,dir1);
101  //cout<<"track1 ";point1.Print();
102  plane1.setO(point1);
103  plane1.setNormal(dir1);
104  rep1->GFAbsTrackRep::extrapolate(plane1);
105  rep2->extrapolateToPoint(point1,point2,norm2);
106  plane2.setO(point1);
107  plane2.setNormal(norm2);
108  rep2->GFAbsTrackRep::extrapolate(plane2);
109 
110  // three points -> fit a parabola -> minimum
111  h=1.;
112  double d0=trkDist(rep1,rep2);
113  //cout<<"dist"<<d0<<endl;
114  rep1->stepalong(h,point1,dir1);
115  //cout<<"track1 ";point1.Print();
116  plane1.setO(point1);
117  plane1.setNormal(dir1);
118  rep1->GFAbsTrackRep::extrapolate(plane1);
119  rep2->extrapolateToPoint(point1,point2,norm2);
120  plane2.setO(point1);
121  plane2.setNormal(norm2);
122  rep2->GFAbsTrackRep::extrapolate(plane2);
123  double d1=trkDist(rep1,rep2);
124  //cout<<"dist"<<d1<<endl;
125  rep1->stepalong(-2.*h,point1,dir1);
126  //cout<<"track1 ";point1.Print();
127  plane1.setO(point1);
128  plane1.setNormal(dir1);
129  rep1->GFAbsTrackRep::extrapolate(plane1);
130  rep2->extrapolateToPoint(point1,point2,norm2);
131  plane2.setO(point1);
132  plane2.setNormal(norm2);
133  rep2->GFAbsTrackRep::extrapolate(plane2);
134  double d2=trkDist(rep1,rep2);
135  //cout<<"dist"<<d2<<endl;
136 
137  double s=(d2-d1)*h;
138  double denom=2.*(d1+d2-2*d0);
139  if(denom==0){
140  TVector3 posA=rep1->GFAbsTrackRep::getPos();
141  TVector3 posB=rep2->GFAbsTrackRep::getPos();
142  TVector3 poca=(posB+posA)*0.5;
143  return poca;
144  // return trkDist(rep1,rep2);
145  }
146  s/=denom;
147 
148  //std:://cout<<"d0="<<d0<<" d1="<<d1<<" d2="<<d2<<" s="<<s<<std::endl;
149 
150  rep1->stepalong(s+h,point1,dir1);
151  //cout<<"track1 ";point1.Print();
152  plane1.setO(point1);
153  plane1.setNormal(dir1);
154  rep1->GFAbsTrackRep::extrapolate(plane1);
155  rep2->extrapolateToPoint(point1,point2,norm2);
156  plane2.setO(point1);
157  plane2.setNormal(norm2);
158  rep2->GFAbsTrackRep::extrapolate(plane2);
159  //double d3=trkDist(rep1,rep2);
160  //cout<<"dist"<<d3<<endl;
161  //cout<<"GFTrackProximity end ########################################################################################################################"<<endl;
162  TVector3 posA=rep1->GFAbsTrackRep::getPos();
163  TVector3 posB=rep2->GFAbsTrackRep::getPos();
164 
165  TVector3 poca=(posB+posA)*0.5;
166  return poca;
167  //return trkDist(rep1,rep2);
168 }
169 
170