00001 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
00002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00003 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00005 #include "FWCore/Utilities/interface/Exception.h"
00006 #include "MagneticField/Engine/interface/MagneticField.h"
00007
00008 using namespace std;
00009
00010 namespace {
00011 inline GlobalPoint mean ( pair<GlobalPoint, GlobalPoint> pr ) {
00012 return GlobalPoint ( (pr.first.x() + pr.second.x() ) / 2. ,
00013 (pr.first.y() + pr.second.y() ) / 2. ,
00014 (pr.first.z() + pr.second.z() ) / 2. );
00015 }
00016
00017 inline double dist ( pair<GlobalPoint, GlobalPoint> pr ) {
00018 return ( pr.first - pr.second ).mag();
00019 }
00020 }
00021
00022 double TwoTrackMinimumDistance::firstAngle() const
00023 {
00024 if (!status_)
00025 throw cms::Exception("TrackingTools/PatternTools","TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
00026 switch ( theCharge ) {
00027 case (hh): return theTTMDhh.firstAngle(); break;
00028 case (hl): return theTTMDhl.firstAngle(); break;
00029 case (ll): return theTTMDll.firstAngle(); break;
00030 }
00031 }
00032
00033 double TwoTrackMinimumDistance::secondAngle() const
00034 {
00035 if (!status_)
00036 throw cms::Exception("TrackingTools/PatternTools","TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
00037 switch ( theCharge ) {
00038 case (hh): return theTTMDhh.secondAngle(); break;
00039 case (hl): return theTTMDhl.secondAngle(); break;
00040 case (ll): return theTTMDll.secondAngle(); break;
00041 }
00042 }
00043
00044
00045 pair <double, double> TwoTrackMinimumDistance::pathLength() const
00046 {
00047 if (!status_)
00048 throw cms::Exception("TrackingTools/PatternTools","TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
00049 switch ( theCharge ) {
00050 case (hh): return theTTMDhh.pathLength(); break;
00051 case (hl): return theTTMDhl.pathLength(); break;
00052 case (ll): return theTTMDll.pathLength(); break;
00053 }
00054 }
00055
00056 pair<GlobalPoint, GlobalPoint> TwoTrackMinimumDistance::points() const
00057 {
00058 if (!status_)
00059 throw cms::Exception("TrackingTools/PatternTools","TwoTrackMinimumDistance::could not compute track crossing. Check status before calling this method!");
00060 return points_;
00061 }
00062
00063 bool
00064 TwoTrackMinimumDistance::calculate(const TrajectoryStateOnSurface & sta,
00065 const TrajectoryStateOnSurface & stb)
00066 {
00067 return calculate ( sta.globalParameters(), stb.globalParameters() );
00068 }
00069
00070
00071 bool
00072 TwoTrackMinimumDistance::calculate(const FreeTrajectoryState & sta,
00073 const FreeTrajectoryState & stb)
00074 {
00075
00076 return calculate ( sta.parameters(), stb.parameters() );
00077 }
00078
00079 bool
00080 TwoTrackMinimumDistance::calculate(const GlobalTrajectoryParameters & sta,
00081 const GlobalTrajectoryParameters & stb)
00082 {
00083 if ((sta.magneticField().inTesla(sta.position()).z() == 0.)||
00084 (stb.magneticField().inTesla(stb.position()).z() == 0.)) {
00085 status_ = pointsLineLine(sta, stb);
00086 } else if ( sta.charge() != 0. && stb.charge() != 0. ) {
00087 status_ = pointsHelixHelix(sta, stb);
00088 } else if ( sta.charge() == 0. && stb.charge() == 0. ) {
00089 status_ = pointsLineLine(sta, stb);
00090 } else {
00091 status_ = pointsHelixLine(sta, stb);
00092 }
00093 return status_;
00094 }
00095
00096 bool
00097 TwoTrackMinimumDistance::pointsLineLine(const GlobalTrajectoryParameters & sta,
00098 const GlobalTrajectoryParameters & stb)
00099 {
00100 theCharge = ll;
00101 if (theTTMDll.calculate(sta, stb)) return false;
00102 points_ = theTTMDll.points();
00103 return true;
00104 }
00105
00106 bool
00107 TwoTrackMinimumDistance::pointsHelixLine(const GlobalTrajectoryParameters & sta,
00108 const GlobalTrajectoryParameters & stb)
00109 {
00110 theCharge = hl;
00111 if (theTTMDhl.calculate(sta, stb, 0.000001)) return false;
00112 points_ = theTTMDhl.points();
00113 return true;
00114 }
00115
00116 bool
00117 TwoTrackMinimumDistance::pointsHelixHelix(const GlobalTrajectoryParameters & sta,
00118 const GlobalTrajectoryParameters & stb)
00119 {
00120 if ( ( sta.position() - stb.position() ).mag() < 1e-7 &&
00121 ( sta.momentum() - stb.momentum() ).mag() < 1e-7 )
00122 {
00123 edm::LogWarning ( "TwoTrackMinimumDistance") << "comparing track with itself!";
00124 };
00125 theCharge = hh;
00126 if ( theModus == FastMode )
00127 {
00128
00129 if ( !(theTTMDhh.calculate ( sta, stb, .0001 )) )
00130 {
00131 points_ = theTTMDhh.points();
00132 return true;
00133 };
00134 };
00135
00136
00137 bool cairStat = theIniAlgo.calculate ( sta, stb );
00138
00139 if (!cairStat) {
00140 edm::LogWarning ( "TwoTrackMinimumDistance" ) << "Computation HelixHelix::CAIR failed.";
00141 if ( theModus == SlowMode ) {
00142 if ( !(theTTMDhh.calculate ( sta, stb, .0001 )) ) {
00143 points_ = theTTMDhh.points();
00144 return true;
00145 }
00146 };
00147
00148 if ( !(theTTMDhh.calculate ( sta, stb, .1 )) ) {
00149 points_ = theTTMDhh.points();
00150 return true;
00151 }
00152 return false;
00153 edm::LogWarning ( "TwoTrackMinimumDistance" ) << "TwoTrackMinimumDistanceHelixHelix failed";
00154 };
00155
00156 pair<GlobalTrajectoryParameters, GlobalTrajectoryParameters >
00157 ini = theIniAlgo.trajectoryParameters();
00158
00159 pair<GlobalPoint, GlobalPoint> inip ( ini.first.position(),
00160 ini.second.position() );
00161 if ( theTTMDhh.calculate ( ini.first, ini.second, .0001 ) ) {
00162 points_ = inip;
00163 } else {
00164 points_ = theTTMDhh.points();
00165
00166 if ( dist ( points_ ) > dist ( inip ) ) points_ = inip;
00167 };
00168 return true;
00169 }
00170
00171
00172 GlobalPoint TwoTrackMinimumDistance::crossingPoint() const
00173 {
00174 return mean ( points_ );
00175 }
00176
00177
00178 float TwoTrackMinimumDistance::distance() const
00179 {
00180 return dist ( points_ );
00181 }