MathTools/src/HelixIntersector.cc

// $Id: HelixIntersector.cc,v 1.3 2009/03/20 13:33:19 loizides Exp $

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

ClassImp(mithep::HelixIntersector)

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.4
Committed:	Mon Jul 20 04:55:44 2009 UTC (15 years, 9 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_032, Mit_031, Mit_025c_branch2, Mit_025c_branch1, Mit_030, Mit_029c, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_025c_branch0, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, V07-05-00, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, ConvRejection-10-06-09, Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, Mit_012g, Mit_012f, Mit_012e, Mit_012d, Mit_012c, Mit_012b, Mit_012a, Mit_012, Mit_011a, Mit_011, Mit_010a, Mit_010, Mit_008, HEAD
Branch point for:	Mit_025c_branch
Changes since 1.3:	+3 -1 lines
Log Message:	Changes for docu
#	User	Rev	Content
1	loizides	1.4	// $Id: HelixIntersector.cc,v 1.3 2009/03/20 13:33:19 loizides Exp $
2	loizides	1.1
3			#include "MitCommon/MathTools/interface/HelixIntersector.h"
4			#include <iostream>
5
6	loizides	1.4	ClassImp(mithep::HelixIntersector)
7
8	loizides	1.1	using namespace std;
9			using namespace mithep;
10
11			//--------------------------------------------------------------------------------------------------
12			HelixIntersector::HelixIntersector(const TVectorD tr1, const TVector3 momentum1,
13			const TVectorD tr2, const TVector3 momentum2) :
14			fISec0(tr1,momentum1,tr2,momentum2),
15			fISec1(tr1,momentum1,tr2,momentum2)
16			{
17	loizides	1.3	// Constructor.
18
19	loizides	1.1	// short cuts
20			TrackParams &trk0 = fISec0.fTrk0;
21			TrackParams &trk1 = fISec0.fTrk1;
22			fISecs[0] = &fISec0;
23			fISecs[1] = &fISec1;
24
25			// Calculate intersections and test for existence
26
27			// line connecting the two centers
28			TVector3 connector = trk1.Center() - trk0.Center();
29
30			// intersection positions relative to connector
31			double para = (connector.Mag2()+trk0.Radius()*trk0.Radius()
32			-trk1.Radius()trk1.Radius())/(2.0connector.Mag());
33			double perpsqr = trk0.Radius()trk0.Radius() - parapara;
34
35			// intersection exist if perp is real -> perp^2 is positive
36			fHasIntersections = (perpsqr >= 0.0);
37
38			// fill Intersection objects appropriately
39			if (fHasIntersections) {
40
41			// vector perpendicular to connector (sign doesn't matter)
42			TVector3 perpvec(connector.y(),-connector.x(),0.0);
43			perpvec.SetMag(sqrt(perpsqr));
44
45			// make both combinations
46			for (int i = 0; i < 2; i++) {
47			// the intersection
48			fISecs[i]->fLocation =
49			trk0.Center() + para * connector.Unit() + double(2i-1) perpvec;
50
51			// calculate tracks/intersection specific quantities for real intersections both tracks have
52			// the same point
53			fISecs[i]->SetLocation(fISecs[i]->fLocation,fISecs[i]->fLocation);
54			}
55
56			if (fabs(fISecs[0]->DeltaZ()) > fabs(fISecs[1]->DeltaZ())) {
57			Intersection* tmp = fISecs[1];
58			fISecs[1] = fISecs[0];
59			fISecs[0] = tmp;
60			}
61
62			}
63			else {
64			// no intersection so we pick points of closest approach
65			TVector3 locTrk0 = trk0.Center() + trk0.Radius() * connector.Unit();
66			TVector3 locTrk1 = trk1.Center() - trk1.Radius() * connector.Unit();
67
68			// mean is the "nominal" location
69			fISecs[0]->fLocation = (locTrk0 + locTrk1) * 0.5;
70
71			// for closest approach points actually on the tracks are used calculate tracks/intersection
72			// specific quantities
73			fISecs[0]->SetLocation(locTrk0,locTrk1);
74			}
75			}
76
77			//--------------------------------------------------------------------------------------------------
78			HelixIntersector::~HelixIntersector()
79			{
80	loizides	1.3	// Destructor.
81	loizides	1.1	}
82
83	loizides	1.3	//--------------------------------------------------------------------------------------------------
84	loizides	1.1	HelixIntersector::TrackParams::TrackParams(const TVectorD trk, const TVector3 momentum) :
85			Helix ((trk)(0),(trk)(1),(trk)(2),(trk)(3),Angle((*trk)(4))),
86			fTrkMom(*momentum)
87			{
88	loizides	1.3	// Calculate center of r-phi plane circles.
89
90	loizides	1.1	double rho = Radius() + Helicity() * D0();
91			double phi = Helicity() * M_PI/2.0 + Phi0();
92			double z = 0.0;
93	loizides	1.3	if (0 && rho < 0) {
94			printf("Cylindrical coordinates supplied with negative Rho\n");
95			printf("Radius = %5f, Helicity = %5f, D0 = %5f, Rho=%5f\n", Radius(), Helicity(), D0(), rho);
96			}
97	loizides	1.1	fCenter.SetXYZ(rhocos(phi),rhosin(phi),z);
98
99			// Original code ....
100			//fCenter.SetRhoPhiZ(Radius() + Helicity() * D0(),
101			// Helicity() * M_PI/2.0 + Phi0(),0.0);
102			}
103
104	loizides	1.3	//--------------------------------------------------------------------------------------------------
105	loizides	1.1	void HelixIntersector::TrackParams::SetLocation(TVector3& location)
106			{
107	loizides	1.3	// Rotation about circle center.
108
109	loizides	1.1	double deltaphi = M_PI - fCenter.Angle(location-fCenter);
110
111			// check for backward locations
112			if (fTrkMom.Dot(location)<0.0)
113			deltaphi = - deltaphi;
114
115
116			// calculate arc length to and z of intersection
117			fArcLen = Radius()*deltaphi/SinTheta();
118			fZAtISec = Z0() + Radius()deltaphiCotTheta();
119
120			// calculate momentum of track at intersection
121			fMomentum = Direction(fArcLen) * fTrkMom.Mag();
122			}
123
124	loizides	1.3	//--------------------------------------------------------------------------------------------------
125	loizides	1.1	HelixIntersector::Intersection::Intersection(const TVectorD tr1, const TVector3 momentum1,
126			const TVectorD tr2, const TVector3 momentum2) :
127			fDeltaZ(0.0),
128			fTrk0 (tr1,momentum1),
129			fTrk1 (tr2,momentum2)
130			{
131	loizides	1.3	// Constructor.
132
133	loizides	1.1	fTrks[0] = &fTrk0;
134			fTrks[1] = &fTrk1;
135			}
136
137	loizides	1.3	//--------------------------------------------------------------------------------------------------
138	loizides	1.1	void HelixIntersector::Intersection::SetLocation(TVector3& loc1, TVector3& loc2)
139			{
140	loizides	1.3	// Calculate tracks specific quantities.
141
142	loizides	1.1	fTrk0.SetLocation(loc1);
143			fTrk1.SetLocation(loc2);
144
145			// combined variables
146			fDeltaZ = fTrk0.fZAtISec - fTrk1.fZAtISec;
147			fLocation.SetZ((fTrk0.fZAtISec + fTrk1.fZAtISec)/2.0);
148			fMomentum = fTrk0.fMomentum + fTrk1.fMomentum;
149			}