PatternTools/test/ClosestApproachInRPhi_t.cpp

0001 #include "TrackingTools/PatternTools/interface/ClosestApproachInRPhi.h"
0002 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistanceHelixHelix.h"
0003
0004 // #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0005 #include "MagneticField/Engine/interface/MagneticField.h"
0006
0007 // #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
0008 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0009
0010 #include "TrackingTools/PatternTools/src/ClosestApproachInRPhi.cc"
0011
0012 #include <iostream>
0013
0014 class ConstMagneticField final : public MagneticField {
0015 public:
0016   ConstMagneticField() { setNominalValue(); }
0017   GlobalVector inTesla(const GlobalPoint&) const override { return GlobalVector(0, 0, 4); }
0018 };
0019
0020 namespace {
0021   inline GlobalPoint mean(std::pair<GlobalPoint, GlobalPoint> pr) {
0022     return GlobalPoint(0.5 * (pr.first.basicVector() + pr.second.basicVector()));
0023   }
0024
0025   inline double dist(std::pair<GlobalPoint, GlobalPoint> pr) { return (pr.first - pr.second).mag(); }
0026 }  // namespace
0027
0028 void compute(GlobalTrajectoryParameters const& gtp1, GlobalTrajectoryParameters const& gtp2) {
0029   ClosestApproachInRPhi ca;
0030   TwoTrackMinimumDistanceHelixHelix TTMDhh;
0031
0032   std::cout << "CAIR" << std::endl;
0033   bool ok = ca.calculate(gtp1, gtp2);
0034   if (!ok)
0035     std::cout << "no intercept!" << std::endl;
0036   else {
0037     std::cout << "distance, xpoint " << ca.distance() << ca.crossingPoint() << std::endl;
0038     std::pair<GlobalTrajectoryParameters, GlobalTrajectoryParameters> tca = ca.trajectoryParameters();
0039     std::cout << tca.first << std::endl;
0040     std::cout << tca.second << std::endl;
0041   }
0042
0043   std::cout << "TTMDhh" << std::endl;
0044   bool nok = TTMDhh.calculate(gtp1, gtp2, .0001);
0045   if (nok)
0046     std::cout << "no intercept!" << std::endl;
0047   else {
0048     std::pair<GlobalPoint, GlobalPoint> pr = TTMDhh.points();
0049     std::cout << "distance, xpoint " << dist(pr) << mean(pr) << std::endl;
0050     // std::pair <GlobalTrajectoryParameters, GlobalTrajectoryParameters > thh = TTMDhh.trajectoryParameters();
0051     // std::cout << thh.first << std::endl;
0052     // std::cout << thh.second << std::endl;
0053   }
0054 }
0055
0056 namespace test {
0057   namespace ClosestApproachInRPhi_t {
0058     int test() {
0059       MagneticField* field = new ConstMagneticField;
0060
0061       {
0062         // going back and forth gtp2 should be identical to gpt1....
0063         GlobalPoint gp1(1, 0, 0);
0064         GlobalVector gv1(1, 1, -1);
0065         GlobalTrajectoryParameters gtp1(gp1, gv1, 1, field);
0066         double bz = field->inTesla(gp1).z() * 2.99792458e-3;
0067         GlobalPoint np(0.504471, -0.498494, 0.497014);
0068         GlobalTrajectoryParameters gtpN = ClosestApproachInRPhi::newTrajectory(np, gtp1, bz);
0069         GlobalTrajectoryParameters gtp2 = ClosestApproachInRPhi::newTrajectory(gp1, gtpN, bz);
0070         std::cout << gtp1 << std::endl;
0071         std::cout << gtpN << std::endl;
0072         std::cout << gtp2 << std::endl;
0073         std::cout << std::endl;
0074       }
0075
0076       {
0077         std::cout << "opposite sign, same direction, same origin: the two circles are tangent to each other at gp1\n"
0078                   << std::endl;
0079         GlobalPoint gp1(0, 0, 0);
0080         GlobalVector gv1(1, 1, 1);
0081         GlobalTrajectoryParameters gtp1(gp1, gv1, 1, field);
0082
0083         GlobalPoint gp2(0, 0, 0);
0084         GlobalVector gv2(1, 1, -1);
0085         GlobalTrajectoryParameters gtp2(gp2, gv2, -1, field);
0086
0087         compute(gtp1, gtp2);
0088         std::cout << std::endl;
0089       }
0090       {
0091         std::cout << " not crossing: the pcas are on the line connecting the two centers\n"
0092                   << "the momenta at the respective pcas shall be parallel as they are perpendicular to the same line\n"
0093                   << "(the one connecting the two centers)\n"
0094                   << std::endl;
0095         GlobalPoint gp1(-1, 0, 0);
0096         GlobalVector gv1(1, 1, 1);
0097         GlobalTrajectoryParameters gtp1(gp1, gv1, -1, field);
0098
0099         GlobalPoint gp2(1, 0, 0);
0100         GlobalVector gv2(1, 1, -1);
0101         GlobalTrajectoryParameters gtp2(gp2, gv2, 1, field);
0102
0103         compute(gtp1, gtp2);
0104         std::cout << std::endl;
0105       }
0106       {
0107         std::cout << "crossing (only opposite changes w.r.t. previous)\n" << std::endl;
0108         GlobalPoint gp1(-1, 0, 0);
0109         GlobalVector gv1(1, 1, 1);
0110         GlobalTrajectoryParameters gtp1(gp1, gv1, 1, field);
0111
0112         GlobalPoint gp2(1, 0, 0);
0113         GlobalVector gv2(1, 1, -1);
0114         GlobalTrajectoryParameters gtp2(gp2, gv2, -1, field);
0115
0116         compute(gtp1, gtp2);
0117         std::cout << std::endl;
0118       }
0119
0120       {
0121         std::cout << "crossing\n" << std::endl;
0122         GlobalPoint gp1(-1, 0, 0);
0123         GlobalVector gv1(1, 1, 1);
0124         GlobalTrajectoryParameters gtp1(gp1, gv1, -1, field);
0125
0126         GlobalPoint gp2(1, 0, 0);
0127         GlobalVector gv2(-1, 1, -1);
0128         GlobalTrajectoryParameters gtp2(gp2, gv2, 1, field);
0129
0130         compute(gtp1, gtp2);
0131         std::cout << std::endl;
0132       }
0133
0134       return 0;
0135     }
0136   }  // namespace ClosestApproachInRPhi_t
0137 }  // namespace test
0138
0139 int main() { return test::ClosestApproachInRPhi_t::test(); }