MathTools/src/HelixIntersector.cc

// $Id: HelixIntersector.cc,v 1.2 2008/11/12 18:22:13 bendavid Exp $

#include "MitCommon/MathTools/interface/HelixIntersector.h"
#include <iostream>

using namespace std;
using namespace mithep;

//--------------------------------------------------------------------------------------------------
HelixIntersector::HelixIntersector(const TVectorD *tr1, const TVector3 *momentum1,
                                   const TVectorD *tr2, const TVector3 *momentum2) :
  fISec0(tr1,momentum1,tr2,momentum2),
  fISec1(tr1,momentum1,tr2,momentum2)
{
  // Constructor.

  // short cuts
  TrackParams &trk0 = fISec0.fTrk0;
  TrackParams &trk1 = fISec0.fTrk1;
  fISecs[0]         = &fISec0;
  fISecs[1]         = &fISec1;

  // Calculate intersections and test for existence

  // line connecting the two centers
  TVector3 connector =  trk1.Center() - trk0.Center();
  
  // intersection positions relative to connector
  double para  = (connector.Mag2()+trk0.Radius()*trk0.Radius()
                  -trk1.Radius()*trk1.Radius())/(2.0*connector.Mag());
  double perpsqr = trk0.Radius()*trk0.Radius() - para*para;
  
  // intersection exist if perp is real -> perp^2 is positive
  fHasIntersections = (perpsqr >= 0.0);
  
  // fill Intersection objects appropriately
  if (fHasIntersections) {

    // vector perpendicular to connector (sign doesn't matter)
    TVector3 perpvec(connector.y(),-connector.x(),0.0);
    perpvec.SetMag(sqrt(perpsqr));
    
    // make both combinations
    for (int i = 0; i < 2; i++) {
      // the intersection
      fISecs[i]->fLocation =
        trk0.Center() + para * connector.Unit() + double(2*i-1) * perpvec;
      
      // calculate tracks/intersection specific quantities for real intersections both tracks have
      // the same point
      fISecs[i]->SetLocation(fISecs[i]->fLocation,fISecs[i]->fLocation);
    }

    if (fabs(fISecs[0]->DeltaZ()) > fabs(fISecs[1]->DeltaZ())) {
      Intersection* tmp = fISecs[1];
      fISecs[1] = fISecs[0];
      fISecs[0] = tmp;
    }

  }
  else {
    // no intersection so we pick points of closest approach
    TVector3 locTrk0 = trk0.Center() + trk0.Radius() * connector.Unit();
    TVector3 locTrk1 = trk1.Center() - trk1.Radius() * connector.Unit();

    // mean is the "nominal" location
    fISecs[0]->fLocation = (locTrk0 + locTrk1) * 0.5;

    // for closest approach points actually on the tracks are used calculate tracks/intersection
    // specific quantities
    fISecs[0]->SetLocation(locTrk0,locTrk1);
  }  
}

//--------------------------------------------------------------------------------------------------
HelixIntersector::~HelixIntersector()
{
  // Destructor.
}

//--------------------------------------------------------------------------------------------------
HelixIntersector::TrackParams::TrackParams(const TVectorD *trk, const TVector3 *momentum) :
  Helix  ((*trk)(0),(*trk)(1),(*trk)(2),(*trk)(3),Angle((*trk)(4))),
  fTrkMom(*momentum)
{
  // Calculate center of r-phi plane circles.

  double rho = Radius() +  Helicity() * D0();
  double phi = Helicity() * M_PI/2.0 + Phi0();
  double z   = 0.0;
  if (0 && rho < 0) {
    printf("Cylindrical coordinates supplied with negative Rho\n");
    printf("Radius = %5f, Helicity = %5f, D0 = %5f, Rho=%5f\n", Radius(), Helicity(), D0(), rho);
  }
  fCenter.SetXYZ(rho*cos(phi),rho*sin(phi),z);

  // Original code ....
  //fCenter.SetRhoPhiZ(Radius() +  Helicity() * D0(),  
  //                   Helicity() * M_PI/2.0 + Phi0(),0.0);
}

//--------------------------------------------------------------------------------------------------
void HelixIntersector::TrackParams::SetLocation(TVector3& location)
{
  // Rotation about circle center.

  double deltaphi = M_PI - fCenter.Angle(location-fCenter);
  
  // check for backward locations
  if (fTrkMom.Dot(location)<0.0)
    deltaphi = - deltaphi;
    

  // calculate arc length to and z of intersection
  fArcLen = Radius()*deltaphi/SinTheta();
  fZAtISec = Z0() + Radius()*deltaphi*CotTheta();

  // calculate momentum of track at intersection
  fMomentum = Direction(fArcLen) * fTrkMom.Mag();
}

//--------------------------------------------------------------------------------------------------
HelixIntersector::Intersection::Intersection(const TVectorD *tr1, const TVector3 *momentum1,
                                             const TVectorD *tr2, const TVector3 *momentum2) :
  fDeltaZ(0.0),
  fTrk0  (tr1,momentum1),
  fTrk1  (tr2,momentum2)
{
  // Constructor.

  fTrks[0] = &fTrk0;
  fTrks[1] = &fTrk1;
}

//--------------------------------------------------------------------------------------------------
void HelixIntersector::Intersection::SetLocation(TVector3& loc1, TVector3& loc2)
{
  // Calculate tracks specific quantities.

  fTrk0.SetLocation(loc1);
  fTrk1.SetLocation(loc2);
  
  // combined variables
  fDeltaZ = fTrk0.fZAtISec - fTrk1.fZAtISec;
  fLocation.SetZ((fTrk0.fZAtISec + fTrk1.fZAtISec)/2.0);
  fMomentum = fTrk0.fMomentum + fTrk1.fMomentum;
}
Revision:	1.3
Committed:	Fri Mar 20 13:33:19 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_009c, Mit_009b, Mit_009a, Mit_009
Changes since 1.2:	+20 -39 lines
Log Message:	Cleanup
#	User	Rev	Content
1	loizides	1.3	// $Id: HelixIntersector.cc,v 1.2 2008/11/12 18:22:13 bendavid Exp $
2	loizides	1.1
3			#include "MitCommon/MathTools/interface/HelixIntersector.h"
4			#include <iostream>
5
6			using namespace std;
7			using namespace mithep;
8
9			//--------------------------------------------------------------------------------------------------
10			HelixIntersector::HelixIntersector(const TVectorD tr1, const TVector3 momentum1,
11			const TVectorD tr2, const TVector3 momentum2) :
12			fISec0(tr1,momentum1,tr2,momentum2),
13			fISec1(tr1,momentum1,tr2,momentum2)
14			{
15	loizides	1.3	// Constructor.
16
17	loizides	1.1	// short cuts
18			TrackParams &trk0 = fISec0.fTrk0;
19			TrackParams &trk1 = fISec0.fTrk1;
20			fISecs[0] = &fISec0;
21			fISecs[1] = &fISec1;
22
23			// Calculate intersections and test for existence
24
25			// line connecting the two centers
26			TVector3 connector = trk1.Center() - trk0.Center();
27
28			// intersection positions relative to connector
29			double para = (connector.Mag2()+trk0.Radius()*trk0.Radius()
30			-trk1.Radius()trk1.Radius())/(2.0connector.Mag());
31			double perpsqr = trk0.Radius()trk0.Radius() - parapara;
32
33			// intersection exist if perp is real -> perp^2 is positive
34			fHasIntersections = (perpsqr >= 0.0);
35
36			// fill Intersection objects appropriately
37			if (fHasIntersections) {
38
39			// vector perpendicular to connector (sign doesn't matter)
40			TVector3 perpvec(connector.y(),-connector.x(),0.0);
41			perpvec.SetMag(sqrt(perpsqr));
42
43			// make both combinations
44			for (int i = 0; i < 2; i++) {
45			// the intersection
46			fISecs[i]->fLocation =
47			trk0.Center() + para * connector.Unit() + double(2i-1) perpvec;
48
49			// calculate tracks/intersection specific quantities for real intersections both tracks have
50			// the same point
51			fISecs[i]->SetLocation(fISecs[i]->fLocation,fISecs[i]->fLocation);
52			}
53
54			if (fabs(fISecs[0]->DeltaZ()) > fabs(fISecs[1]->DeltaZ())) {
55			Intersection* tmp = fISecs[1];
56			fISecs[1] = fISecs[0];
57			fISecs[0] = tmp;
58			}
59
60			}
61			else {
62			// no intersection so we pick points of closest approach
63			TVector3 locTrk0 = trk0.Center() + trk0.Radius() * connector.Unit();
64			TVector3 locTrk1 = trk1.Center() - trk1.Radius() * connector.Unit();
65
66			// mean is the "nominal" location
67			fISecs[0]->fLocation = (locTrk0 + locTrk1) * 0.5;
68
69			// for closest approach points actually on the tracks are used calculate tracks/intersection
70			// specific quantities
71			fISecs[0]->SetLocation(locTrk0,locTrk1);
72			}
73			}
74
75			//--------------------------------------------------------------------------------------------------
76			HelixIntersector::~HelixIntersector()
77			{
78	loizides	1.3	// Destructor.
79	loizides	1.1	}
80
81	loizides	1.3	//--------------------------------------------------------------------------------------------------
82	loizides	1.1	HelixIntersector::TrackParams::TrackParams(const TVectorD trk, const TVector3 momentum) :
83			Helix ((trk)(0),(trk)(1),(trk)(2),(trk)(3),Angle((*trk)(4))),
84			fTrkMom(*momentum)
85			{
86	loizides	1.3	// Calculate center of r-phi plane circles.
87
88	loizides	1.1	double rho = Radius() + Helicity() * D0();
89			double phi = Helicity() * M_PI/2.0 + Phi0();
90			double z = 0.0;
91	loizides	1.3	if (0 && rho < 0) {
92			printf("Cylindrical coordinates supplied with negative Rho\n");
93			printf("Radius = %5f, Helicity = %5f, D0 = %5f, Rho=%5f\n", Radius(), Helicity(), D0(), rho);
94			}
95	loizides	1.1	fCenter.SetXYZ(rhocos(phi),rhosin(phi),z);
96
97			// Original code ....
98			//fCenter.SetRhoPhiZ(Radius() + Helicity() * D0(),
99			// Helicity() * M_PI/2.0 + Phi0(),0.0);
100			}
101
102	loizides	1.3	//--------------------------------------------------------------------------------------------------
103	loizides	1.1	void HelixIntersector::TrackParams::SetLocation(TVector3& location)
104			{
105	loizides	1.3	// Rotation about circle center.
106
107	loizides	1.1	double deltaphi = M_PI - fCenter.Angle(location-fCenter);
108
109			// check for backward locations
110			if (fTrkMom.Dot(location)<0.0)
111			deltaphi = - deltaphi;
112
113
114			// calculate arc length to and z of intersection
115			fArcLen = Radius()*deltaphi/SinTheta();
116			fZAtISec = Z0() + Radius()deltaphiCotTheta();
117
118			// calculate momentum of track at intersection
119			fMomentum = Direction(fArcLen) * fTrkMom.Mag();
120			}
121
122	loizides	1.3	//--------------------------------------------------------------------------------------------------
123	loizides	1.1	HelixIntersector::Intersection::Intersection(const TVectorD tr1, const TVector3 momentum1,
124			const TVectorD tr2, const TVector3 momentum2) :
125			fDeltaZ(0.0),
126			fTrk0 (tr1,momentum1),
127			fTrk1 (tr2,momentum2)
128			{
129	loizides	1.3	// Constructor.
130
131	loizides	1.1	fTrks[0] = &fTrk0;
132			fTrks[1] = &fTrk1;
133			}
134
135	loizides	1.3	//--------------------------------------------------------------------------------------------------
136	loizides	1.1	void HelixIntersector::Intersection::SetLocation(TVector3& loc1, TVector3& loc2)
137			{
138	loizides	1.3	// Calculate tracks specific quantities.
139
140	loizides	1.1	fTrk0.SetLocation(loc1);
141			fTrk1.SetLocation(loc2);
142
143			// combined variables
144			fDeltaZ = fTrk0.fZAtISec - fTrk1.fZAtISec;
145			fLocation.SetZ((fTrk0.fZAtISec + fTrk1.fZAtISec)/2.0);
146			fMomentum = fTrk0.fMomentum + fTrk1.fMomentum;
147			}