RungeKutta/test/RKTest.cc

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "MagneticField/Engine/interface/MagneticField.h"
0010 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0011
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014
0015 #include <iostream>
0016 #include <string>
0017 #include <map>
0018
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020
0021 ////////////////////////////////////////////////////////////////////////////
0022
0023 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0024 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0025 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0026 #include "MagneticField/VolumeGeometry/interface/MagVolumeOutsideValidity.h"
0027 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0028
0029 ////////////////////////////////////////////////////////////////////////////
0030
0031 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0032 #include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
0033
0034 class RKTestField final : public MagneticField {
0035 public:
0036   RKTestField() { setNominalValue(); }
0037   GlobalVector inTesla(const GlobalPoint&) const override { return GlobalVector(0, 0, 4); }
0038 };
0039
0040 using namespace std;
0041
0042 class RKTest : public edm::one::EDAnalyzer<> {
0043 public:
0044   RKTest(const edm::ParameterSet& pset) : magFieldToken(esConsumes()) {}
0045
0046   ~RKTest() = default;
0047
0048   virtual void analyze(const edm::Event& event, const edm::EventSetup& setup) {
0049     using namespace edm;
0050     const ESHandle<MagneticField> magfield = setup.getHandle(magFieldToken);
0051     propagateInCentralVolume(&(*magfield));
0052   }
0053
0054 private:
0055   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken;
0056   typedef TrajectoryStateOnSurface TSOS;
0057
0058   void propagateInCentralVolume(const MagneticField* field) const;
0059   Surface::RotationType rotation(const GlobalVector& zAxis) const;
0060 };
0061
0062 void RKTest::propagateInCentralVolume(const MagneticField* field) const {
0063   // RKTestField is used internally in RKTestPropagator
0064   // In the following code "field" or "&TestField" are interchangeable
0065   // They should give identical field values if "field" is produced by VolumeBasedMagneticField
0066   // with the non-default option "useParametrizedTrackerField = true"
0067
0068   RKTestField TestField;  // Not needed if you want to use field instead of &TestField
0069
0070   //  RKTestPropagator RKprop ( &TestField, alongMomentum );
0071   //AnalyticalPropagator ANprop  ( &TestField, alongMomentum);
0072   defaultRKPropagator::Product prod(field, alongMomentum, 5.e-5);
0073   auto& RKprop = prod.propagator;
0074   AnalyticalPropagator ANprop(field, alongMomentum);
0075
0076   for (float phi = -3.14; phi < 3.14; phi += 0.5) {
0077     for (float costh = -0.99; costh < +0.99; costh += 0.3) {
0078       //Define starting position and momentum
0079       float sinth = sqrt(1 - costh * costh);
0080       for (float p = 0.5; p < 12; p += 1.) {
0081         GlobalVector startingMomentum(p * sin(phi) * sinth, p * cos(phi) * sinth, p * costh);
0082         for (float z = -100.; z < 100.; z += 10.) {
0083           GlobalPoint startingPosition(0, 0, z);
0084           //Define starting plane
0085           PlaneBuilder pb;
0086           Surface::RotationType rot = rotation(startingMomentum);
0087           PlaneBuilder::ReturnType startingPlane = pb.plane(startingPosition, rot);
0088           // Define end plane
0089           for (float d = 10.; d < 150.; d += 10.) {
0090             cout << "And now trying costh, phi, p, z, dist = " << costh << ", " << phi << ", " << p << ", " << z << ", "
0091                  << d << endl;
0092
0093             float propDistance = d;  // 100 cm
0094             GlobalPoint targetPos = startingPosition + propDistance * startingMomentum.unit();
0095             PlaneBuilder::ReturnType EndPlane = pb.plane(targetPos, rot);
0096             // Define error matrix
0097             ROOT::Math::SMatrixIdentity id;
0098             AlgebraicSymMatrix55 C(id);
0099             C *= 0.01;
0100             CurvilinearTrajectoryError err(C);
0101
0102             TSOS startingStateP(
0103                 GlobalTrajectoryParameters(startingPosition, startingMomentum, 1, &TestField), err, *startingPlane);
0104
0105             TSOS startingStateM(
0106                 GlobalTrajectoryParameters(startingPosition, startingMomentum, -1, &TestField), err, *startingPlane);
0107
0108             try {
0109               TSOS trackStateP = RKprop.propagate(startingStateP, *EndPlane);
0110               if (trackStateP.isValid())
0111                 cout << "Succesfully finished Positive track propagation  -------------- with RK: "
0112                      << trackStateP.globalPosition() << endl;
0113               TSOS trackStateP2 = ANprop.propagate(startingStateP, *EndPlane);
0114               if (trackStateP2.isValid())
0115                 cout << "Succesfully finished Positive track propagation  -------------- with AN: "
0116                      << trackStateP2.globalPosition() << endl;
0117             } catch (MagVolumeOutsideValidity& duh) {
0118               cout << "MagVolumeOutsideValidity not properly caught!! Lost this track " << endl;
0119             }
0120
0121             try {
0122               TSOS trackStateM = RKprop.propagate(startingStateM, *EndPlane);
0123               if (trackStateM.isValid())
0124                 cout << "Succesfully finished Negative track propagation  -------------- with RK: "
0125                      << trackStateM.globalPosition() << endl;
0126               TSOS trackStateM2 = ANprop.propagate(startingStateM, *EndPlane);
0127               if (trackStateM2.isValid())
0128                 cout << "Succesfully finished Negative track propagation  -------------- with AN: "
0129                      << trackStateM2.globalPosition() << endl;
0130             } catch (MagVolumeOutsideValidity& duh) {
0131               cout << "MagVolumeOutsideValidity not properly caught!! Lost this track " << endl;
0132             }
0133           }
0134         }
0135       }
0136     }
0137   }
0138   cout << " Succesfully reached the END of this test !!!!!!!!!! " << endl;
0139 }
0140
0141 Surface::RotationType RKTest::rotation(const GlobalVector& zDir) const {
0142   GlobalVector zAxis = zDir.unit();
0143   GlobalVector yAxis(zAxis.y(), -zAxis.x(), 0);
0144   GlobalVector xAxis = yAxis.cross(zAxis);
0145   return Surface::RotationType(xAxis, yAxis, zAxis);
0146 }
0147
0148 DEFINE_FWK_MODULE(RKTest);