Geant4e/test/SimpleGeant4ePropagatorTest.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"  //For define_fwk_module
0006
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Utilities/interface/InputTag.h"
0010
0011 //- Timing
0012 //#include "Utilities/Timing/interface/TimingReport.h"
0013
0014 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0015
0016 //- Magnetic field
0017 #include "MagneticField/Engine/interface/MagneticField.h"
0018 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0019 #include "SimG4Core/MagneticField/interface/Field.h"
0020 #include "SimG4Core/MagneticField/interface/FieldBuilder.h"
0021
0022 //- Propagator
0023 #include "TrackPropagation/Geant4e/interface/ConvertFromToCLHEP.h"
0024 #include "TrackPropagation/Geant4e/interface/Geant4ePropagator.h"
0025 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0026 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0027
0028 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0029 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0030 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0031 #include "MagneticField/VolumeGeometry/interface/MagVolumeOutsideValidity.h"
0032 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0033
0034 //- Geant4
0035 #include "G4TransportationManager.hh"
0036
0037 class SimpleGeant4ePropagatorTest final : public edm::one::EDAnalyzer<> {
0038 public:
0039   explicit SimpleGeant4ePropagatorTest(const edm::ParameterSet &);
0040   ~SimpleGeant4ePropagatorTest() override = default;
0041
0042   void analyze(const edm::Event &, const edm::EventSetup &) override;
0043
0044 protected:
0045   // tokens
0046   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken;
0047
0048   Propagator *thePropagator;
0049 };
0050
0051 SimpleGeant4ePropagatorTest::SimpleGeant4ePropagatorTest(const edm::ParameterSet &iConfig)
0052     : magFieldToken(esConsumes()), thePropagator(nullptr) {}
0053
0054 namespace {
0055   Surface::RotationType rotation(const GlobalVector &zDir) {
0056     GlobalVector zAxis = zDir.unit();
0057     GlobalVector yAxis(zAxis.y(), -zAxis.x(), 0);
0058     GlobalVector xAxis = yAxis.cross(zAxis);
0059     return Surface::RotationType(xAxis, yAxis, zAxis);
0060   }
0061 }  // namespace
0062
0063 void SimpleGeant4ePropagatorTest::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0064   using namespace edm;
0065
0066   std::cout << "Starting G4e test..." << std::endl;
0067
0068   ///////////////////////////////////////
0069   // Construct Magnetic Field
0070   const ESHandle<MagneticField> bField = iSetup.getHandle(magFieldToken);
0071   if (bField.isValid())
0072     std::cout << "G4e -- Magnetic field is valid. Value in (0,0,0): " << bField->inTesla(GlobalPoint(0, 0, 0)).mag()
0073               << " Tesla " << std::endl;
0074   else
0075     LogError("Geant4e") << "G4e -- NO valid Magnetic field" << std::endl;
0076
0077   // Initialise the propagator
0078   if (!thePropagator)
0079     thePropagator = new Geant4ePropagator(bField.product());
0080
0081   if (thePropagator)
0082     std::cout << "Propagator built!" << std::endl;
0083   else
0084     LogError("Geant4e") << "Could not build propagator!" << std::endl;
0085
0086   GlobalVector p3T(10., 10., 2.);
0087   std::cout << "*** Phi (rad): " << p3T.phi() << " - Phi(deg)" << p3T.phi().degrees();
0088   std::cout << "Track P.: " << p3T << "\nTrack P.: PT=" << p3T.perp() << "\tEta=" << p3T.eta()
0089             << "\tPhi=" << p3T.phi().degrees() << "--> Rad: Phi=" << p3T.phi() << std::endl;
0090
0091   GlobalPoint r3T(0., 0., 0.);
0092   std::cout << "Init point: " << r3T << "\nInit point Ro=" << r3T.perp() << "\tEta=" << r3T.eta()
0093             << "\tPhi=" << r3T.phi().degrees() << std::endl;
0094
0095   //- Charge
0096   int charge = 1;
0097   std::cout << "Track charge = " << charge << std::endl;
0098
0099   //- Initial covariance matrix is unity 10-6
0100   ROOT::Math::SMatrixIdentity id;
0101   AlgebraicSymMatrix55 C(id);
0102   C *= 0.01;
0103   CurvilinearTrajectoryError covT(C);
0104
0105   PlaneBuilder pb;
0106   Surface::RotationType rot = rotation(p3T);
0107   // Define end planes
0108   for (float d = 50.; d < 700.; d += 50.) {
0109     float propDistance = d;  // 100 cm
0110     std::cout << "G4e -- Extrapolatation distance " << d << " cm" << std::endl;
0111     GlobalPoint targetPos = r3T + propDistance * p3T.unit();
0112     auto endPlane = pb.plane(targetPos, rot);
0113
0114     //- Build FreeTrajectoryState
0115     GlobalTrajectoryParameters trackPars(r3T, p3T, charge, &*bField);
0116     FreeTrajectoryState ftsTrack(trackPars, covT);
0117
0118     // Propagate: Need to explicetely
0119     TrajectoryStateOnSurface tSOSDest = thePropagator->propagate(ftsTrack, *endPlane);
0120     if (!tSOSDest.isValid()) {
0121       std::cout << "TSOS not valid? Propagation failed...." << std::endl;
0122       continue;
0123     }
0124
0125     auto posExtrap = tSOSDest.freeState()->position();
0126     auto momExtrap = tSOSDest.freeState()->momentum();
0127     std::cout << "G4e -- Extrapolated position:" << posExtrap << " cm\n"
0128               << "G4e --       (Rho, eta, phi): (" << posExtrap.perp() << " cm, " << posExtrap.eta() << ", "
0129               << posExtrap.phi() << ')' << std::endl;
0130     std::cout << "G4e -- Extrapolated momentum:" << momExtrap << " GeV\n"
0131               << "G4e --       (pt, eta, phi): (" << momExtrap.perp() << " cm, " << momExtrap.eta() << ", "
0132               << momExtrap.phi() << ')' << std::endl;
0133   }
0134 }
0135
0136 // define this as a plug-in
0137 DEFINE_FWK_MODULE(SimpleGeant4ePropagatorTest);