Geant4e/src/Geant4ePropagator.cc

0001 #include <sstream>
0002
0003 // Geant4e
0004 #include "TrackPropagation/Geant4e/interface/ConvertFromToCLHEP.h"
0005 #include "TrackPropagation/Geant4e/interface/Geant4ePropagator.h"
0006
0007 // CMSSW
0008 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0009 #include "MagneticField/Engine/interface/MagneticField.h"
0010 #include "TrackingTools/TrajectoryState/interface/SurfaceSideDefinition.h"
0011 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0012
0013 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0014 #include "DataFormats/GeometrySurface/interface/Plane.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0017
0018 // Geant4
0019 #include "G4Box.hh"
0020 #include "G4ErrorCylSurfaceTarget.hh"
0021 #include "G4ErrorFreeTrajState.hh"
0022 #include "G4ErrorPlaneSurfaceTarget.hh"
0023 #include "G4ErrorPropagatorData.hh"
0024 #include "G4ErrorRunManagerHelper.hh"
0025 #include "G4EventManager.hh"
0026 #include "G4Field.hh"
0027 #include "G4FieldManager.hh"
0028 #include "G4GeometryTolerance.hh"
0029 #include "G4SteppingControl.hh"
0030 #include "G4TransportationManager.hh"
0031 #include "G4Tubs.hh"
0032 #include "G4UImanager.hh"
0033 #include "G4ErrorPropagationNavigator.hh"
0034 #include "G4RunManagerKernel.hh"
0035 #include "G4StateManager.hh"
0036
0037 // CLHEP
0038 #include <CLHEP/Units/SystemOfUnits.h>
0039
0040 /** Constructor.
0041  */
0042 Geant4ePropagator::Geant4ePropagator(const MagneticField *field,
0043                                      std::string particleName,
0044                                      PropagationDirection dir,
0045                                      double plimit)
0046     : Propagator(dir),
0047       theField(field),
0048       theParticleName(particleName),
0049       theG4eManager(G4ErrorPropagatorManager::GetErrorPropagatorManager()),
0050       theG4eData(G4ErrorPropagatorData::GetErrorPropagatorData()),
0051       plimit_(plimit) {
0052   LogDebug("Geant4e") << "Geant4e Propagator initialized";
0053
0054   // has to be called here, doing it later will not load the G4 physics list
0055   // properly when using the G4 ES Producer. Reason: unclear
0056   ensureGeant4eIsInitilized(true);
0057 }
0058
0059 /** Destructor.
0060  */
0061 Geant4ePropagator::~Geant4ePropagator() {
0062   LogDebug("Geant4e") << "Geant4ePropagator::~Geant4ePropagator()" << std::endl;
0063
0064   // don't close the g4 Geometry here, because the propagator might have been
0065   // cloned
0066   // but there is only one, globally-shared Geometry
0067 }
0068
0069 //
0070 ////////////////////////////////////////////////////////////////////////////
0071 //
0072
0073 /** Propagate from a free state (e.g. position and momentum in
0074  *  in global cartesian coordinates) to a plane.
0075  */
0076
0077 void Geant4ePropagator::ensureGeant4eIsInitilized(bool) const {
0078   LogDebug("Geant4ePropagator") << "G4 propagator starts isInitialized, theField: " << theField;
0079
0080   auto man = G4RunManagerKernel::GetRunManagerKernel();
0081   if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_PreInit) {
0082     man->SetVerboseLevel(0);
0083     theG4eManager->InitGeant4e();
0084
0085     // define 10 mm step limit for propagator
0086     G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/stepLength 10.0 mm");
0087   }
0088   const G4Field *field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
0089   if (field == nullptr) {
0090     edm::LogError("Geant4e") << "No G4 magnetic field defined";
0091   }
0092   LogDebug("Geant4ePropagator") << "G4 propagator initialized; field: " << field;
0093 }
0094
0095 template <>
0096 Geant4ePropagator::ErrorTargetPair Geant4ePropagator::transformToG4SurfaceTarget(const Plane &pDest,
0097                                                                                  bool moveTargetToEndOfSurface) const {
0098   //* Get position and normal (orientation) of the destination plane
0099   GlobalPoint posPlane = pDest.toGlobal(LocalPoint(0, 0, 0));
0100   GlobalVector normalPlane = pDest.toGlobal(LocalVector(0, 0, 1.));
0101   normalPlane = normalPlane.unit();
0102
0103   //* Transform this into HepGeom::Point3D<double>  and
0104   // HepGeom::Normal3D<double>  that define a plane for
0105   //  Geant4e.
0106   //  CMS uses cm and GeV while Geant4 uses mm and MeV
0107   HepGeom::Point3D<double> surfPos = TrackPropagation::globalPointToHepPoint3D(posPlane);
0108   HepGeom::Normal3D<double> surfNorm = TrackPropagation::globalVectorToHepNormal3D(normalPlane);
0109
0110   //* Set the target surface
0111   return ErrorTargetPair(false, std::make_shared<G4ErrorPlaneSurfaceTarget>(surfNorm, surfPos));
0112 }
0113
0114 template <>
0115 Geant4ePropagator::ErrorTargetPair Geant4ePropagator::transformToG4SurfaceTarget(const Cylinder &pDest,
0116                                                                                  bool moveTargetToEndOfSurface) const {
0117   // Get Cylinder parameters.
0118   // CMS uses cm and GeV while Geant4 uses mm and MeV.
0119   // - Radius
0120   G4float radCyl = pDest.radius() * CLHEP::cm;
0121   // - Position: PositionType & GlobalPoint are Basic3DPoint<float,GlobalTag>
0122   G4ThreeVector posCyl = TrackPropagation::globalPointToHep3Vector(pDest.position());
0123   // - Rotation: Type in CMSSW is RotationType == TkRotation<T>, T=float
0124   G4RotationMatrix rotCyl = TrackPropagation::tkRotationFToHepRotation(pDest.rotation());
0125
0126   // DEBUG
0127   TkRotation<float> rotation = pDest.rotation();
0128   LogDebug("Geant4e") << "G4e -  TkRotation" << rotation;
0129   LogDebug("Geant4e") << "G4e -  G4Rotation" << rotCyl << "mm";
0130
0131   return ErrorTargetPair(!moveTargetToEndOfSurface, std::make_shared<G4ErrorCylSurfaceTarget>(radCyl, posCyl, rotCyl));
0132 }
0133
0134 template <>
0135 std::string Geant4ePropagator::getSurfaceType(Cylinder const &c) const {
0136   return "Cylinder";
0137 }
0138
0139 template <>
0140 std::string Geant4ePropagator::getSurfaceType(Plane const &c) const {
0141   return "Plane";
0142 }
0143
0144 std::string Geant4ePropagator::generateParticleName(int charge) const {
0145   std::string particleName = theParticleName;
0146
0147   if (charge > 0) {
0148     particleName += "+";
0149   }
0150   if (charge < 0) {
0151     particleName += "-";
0152   }
0153
0154   LogDebug("Geant4e") << "G4e -  Particle name: " << particleName;
0155
0156   return particleName;
0157 }
0158
0159 template <>
0160 bool Geant4ePropagator::configureAnyPropagation(G4ErrorMode &mode,
0161                                                 Plane const &pDest,
0162                                                 GlobalPoint const &cmsInitPos,
0163                                                 GlobalVector const &cmsInitMom) const {
0164   if (cmsInitMom.mag() < plimit_)
0165     return false;
0166   if (pDest.localZ(cmsInitPos) * pDest.localZ(cmsInitMom) < 0) {
0167     mode = G4ErrorMode_PropForwards;
0168     LogDebug("Geant4e") << "G4e -  Propagator mode is \'forwards\' indirect "
0169                            "via the Any direction"
0170                         << std::endl;
0171   } else {
0172     mode = G4ErrorMode_PropBackwards;
0173     LogDebug("Geant4e") << "G4e -  Propagator mode is \'backwards\' indirect "
0174                            "via the Any direction"
0175                         << std::endl;
0176   }
0177
0178   return true;
0179 }
0180
0181 template <>
0182 bool Geant4ePropagator::configureAnyPropagation(G4ErrorMode &mode,
0183                                                 Cylinder const &pDest,
0184                                                 GlobalPoint const &cmsInitPos,
0185                                                 GlobalVector const &cmsInitMom) const {
0186   if (cmsInitMom.mag() < plimit_)
0187     return false;
0188   //------------------------------------
0189   // For cylinder assume outside is backwards, inside is along
0190   // General use for particles from collisions
0191   LocalPoint lpos = pDest.toLocal(cmsInitPos);
0192   Surface::Side theSide = pDest.side(lpos, 0);
0193   if (theSide == SurfaceOrientation::positiveSide) {  // outside cylinder
0194     mode = G4ErrorMode_PropBackwards;
0195     LogDebug("Geant4e") << "G4e -  Propagator mode is \'backwards\' indirect "
0196                            "via the Any direction";
0197   } else {  // inside cylinder
0198     mode = G4ErrorMode_PropForwards;
0199     LogDebug("Geant4e") << "G4e -  Propagator mode is \'forwards\' indirect "
0200                            "via the Any direction";
0201   }
0202
0203   return true;
0204 }
0205
0206 template <class SurfaceType>
0207 bool Geant4ePropagator::configurePropagation(G4ErrorMode &mode,
0208                                              SurfaceType const &pDest,
0209                                              GlobalPoint const &cmsInitPos,
0210                                              GlobalVector const &cmsInitMom) const {
0211   if (cmsInitMom.mag() < plimit_)
0212     return false;
0213   if (propagationDirection() == oppositeToMomentum) {
0214     mode = G4ErrorMode_PropBackwards;
0215     LogDebug("Geant4e") << "G4e -  Propagator mode is \'backwards\' " << std::endl;
0216   } else if (propagationDirection() == alongMomentum) {
0217     mode = G4ErrorMode_PropForwards;
0218     LogDebug("Geant4e") << "G4e -  Propagator mode is \'forwards\'" << std::endl;
0219   } else if (propagationDirection() == anyDirection) {
0220     if (configureAnyPropagation(mode, pDest, cmsInitPos, cmsInitMom) == false)
0221       return false;
0222   } else {
0223     edm::LogError("Geant4e") << "G4e - Unsupported propagation mode";
0224     return false;
0225   }
0226   return true;
0227 }
0228
0229 template <class SurfaceType>
0230 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateGeneric(const FreeTrajectoryState &ftsStart,
0231                                                                                 const SurfaceType &pDest) const {
0232   ///////////////////////////////
0233   // Construct the target surface
0234   //
0235   //* Set the target surface
0236
0237   ErrorTargetPair g4eTarget_center = transformToG4SurfaceTarget(pDest, false);
0238
0239   // * Get the starting point and direction and convert them to
0240   // CLHEP::Hep3Vector
0241   //   for G4. CMS uses cm and GeV while Geant4 uses mm and MeV
0242   GlobalPoint cmsInitPos = ftsStart.position();
0243   GlobalVector cmsInitMom = ftsStart.momentum();
0244   bool flipped = false;
0245   if (propagationDirection() == oppositeToMomentum) {
0246     // flip the momentum vector as Geant4 will not do this
0247     // on it's own in a backward propagation
0248     cmsInitMom = -cmsInitMom;
0249     flipped = true;
0250   }
0251
0252   // Set the mode of propagation according to the propagation direction
0253   G4ErrorMode mode = G4ErrorMode_PropForwards;
0254   if (!configurePropagation(mode, pDest, cmsInitPos, cmsInitMom))
0255     return TsosPP(TrajectoryStateOnSurface(), 0.0f);
0256
0257   // re-check propagation direction chosen in case of AnyDirection
0258   if (mode == G4ErrorMode_PropBackwards && !flipped)
0259     cmsInitMom = -cmsInitMom;
0260
0261   CLHEP::Hep3Vector g4InitPos = TrackPropagation::globalPointToHep3Vector(cmsInitPos);
0262   CLHEP::Hep3Vector g4InitMom = TrackPropagation::globalVectorToHep3Vector(cmsInitMom * CLHEP::GeV);
0263
0264   debugReportTrackState("intitial", cmsInitPos, g4InitPos, cmsInitMom, g4InitMom, pDest);
0265
0266   // Set the mode of propagation according to the propagation direction
0267   // G4ErrorMode mode = G4ErrorMode_PropForwards;
0268
0269   // if (!configurePropagation(mode, pDest, cmsInitPos, cmsInitMom))
0270   //    return TsosPP(TrajectoryStateOnSurface(), 0.0f);
0271
0272   ///////////////////////////////
0273   // Set the error and trajectories, and finally propagate
0274   //
0275   G4ErrorTrajErr g4error(5, 1);
0276   if (ftsStart.hasError()) {
0277     CurvilinearTrajectoryError initErr;
0278     initErr = ftsStart.curvilinearError();
0279     g4error = TrackPropagation::algebraicSymMatrix55ToG4ErrorTrajErr(initErr, ftsStart.charge());
0280     LogDebug("Geant4e") << "CMS -  Error matrix: " << std::endl << initErr.matrix();
0281   } else {
0282     LogDebug("Geant4e") << "No error matrix available" << std::endl;
0283     return TsosPP(TrajectoryStateOnSurface(), 0.0f);
0284   }
0285
0286   LogDebug("Geant4e") << "G4e -  Error matrix: " << std::endl << g4error;
0287
0288   // in CMSSW, the state errors are deflated when performing the backward
0289   // propagation
0290   if (mode == G4ErrorMode_PropForwards) {
0291     G4ErrorPropagatorData::GetErrorPropagatorData()->SetStage(G4ErrorStage_Inflation);
0292   } else if (mode == G4ErrorMode_PropBackwards) {
0293     G4ErrorPropagatorData::GetErrorPropagatorData()->SetStage(G4ErrorStage_Deflation);
0294   }
0295
0296   G4ErrorFreeTrajState g4eTrajState(generateParticleName(ftsStart.charge()), g4InitPos, g4InitMom, g4error);
0297   LogDebug("Geant4e") << "G4e -  Traj. State: " << (g4eTrajState);
0298
0299   //////////////////////////////
0300   // Propagate
0301   int iterations = 0;
0302   double finalPathLength = 0;
0303
0304   HepGeom::Point3D<double> finalRecoPos;
0305
0306   G4ErrorPropagatorData::GetErrorPropagatorData()->SetMode(mode);
0307
0308   theG4eData->SetTarget(g4eTarget_center.second.get());
0309   LogDebug("Geant4e") << "Running Propagation to the RECO surface" << std::endl;
0310
0311   theG4eManager->InitTrackPropagation();
0312
0313   // re-initialize navigator to avoid mismatches and/or segfaults
0314   theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointAndSetup(
0315       g4InitPos, &g4InitMom, /*pRelativeSearch = */ false, /*ignoreDirection = */ false);
0316
0317   bool continuePropagation = true;
0318   while (continuePropagation) {
0319     iterations++;
0320     LogDebug("Geant4e") << std::endl << "step count " << iterations << " step length " << finalPathLength;
0321
0322     // re-initialize navigator to avoid mismatches and/or segfaults
0323     theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointWithinVolume(g4eTrajState.GetPosition());
0324
0325     const int ierr = theG4eManager->PropagateOneStep(&g4eTrajState, mode);
0326
0327     if (ierr != 0) {
0328       // propagation failed, return invalid track state
0329       return TsosPP(TrajectoryStateOnSurface(), 0.0f);
0330     }
0331
0332     const float thisPathLength = TrackPropagation::g4doubleToCmsDouble(g4eTrajState.GetG4Track()->GetStepLength());
0333
0334     LogDebug("Geant4e") << "step Length was " << thisPathLength << " cm, current global position: "
0335                         << TrackPropagation::hepPoint3DToGlobalPoint(g4eTrajState.GetPosition()) << std::endl;
0336
0337     finalPathLength += thisPathLength;
0338
0339     // if (std::fabs(finalPathLength) > 10000.0f)
0340     if (std::fabs(finalPathLength) > 200.0f) {
0341       LogDebug("Geant4e") << "ERROR: Quitting propagation: path length mega large" << std::endl;
0342       theG4eManager->GetPropagator()->InvokePostUserTrackingAction(g4eTrajState.GetG4Track());
0343       continuePropagation = false;
0344       LogDebug("Geant4e") << "WARNING: Quitting propagation: max path length "
0345                              "exceeded, returning invalid state"
0346                           << std::endl;
0347
0348       // reached maximum path length, bail out
0349       return TsosPP(TrajectoryStateOnSurface(), 0.0f);
0350     }
0351
0352     if (theG4eManager->GetPropagator()->CheckIfLastStep(g4eTrajState.GetG4Track())) {
0353       theG4eManager->GetPropagator()->InvokePostUserTrackingAction(g4eTrajState.GetG4Track());
0354       continuePropagation = false;
0355     }
0356   }
0357
0358   // CMSSW Tracking convention, backward propagations have negative path length
0359   if (propagationDirection() == oppositeToMomentum)
0360     finalPathLength = -finalPathLength;
0361
0362   // store the correct location for the hit on the RECO surface
0363   LogDebug("Geant4e") << "Position on the RECO surface" << g4eTrajState.GetPosition() << std::endl;
0364   finalRecoPos = g4eTrajState.GetPosition();
0365
0366   theG4eManager->EventTermination();
0367
0368   LogDebug("Geant4e") << "Final position of the Track :" << g4eTrajState.GetPosition() << std::endl;
0369
0370   //////////////////////////////
0371   // Retrieve the state in the end from Geant4e, convert them to CMS vectors
0372   // and points, and build global trajectory parameters.
0373   // CMS uses cm and GeV while Geant4 uses mm and MeV
0374   //
0375   const HepGeom::Vector3D<double> momEnd = g4eTrajState.GetMomentum();
0376
0377   // use the hit on the the RECO plane as the final position to be d'accor with
0378   // the RecHit measurements
0379   const GlobalPoint posEndGV = TrackPropagation::hepPoint3DToGlobalPoint(finalRecoPos);
0380   GlobalVector momEndGV = TrackPropagation::hep3VectorToGlobalVector(momEnd) / CLHEP::GeV;
0381
0382   debugReportTrackState("final", posEndGV, finalRecoPos, momEndGV, momEnd, pDest);
0383
0384   // Get the error covariance matrix from Geant4e. It comes in curvilinear
0385   // coordinates so use the appropiate CMS class
0386   G4ErrorTrajErr g4errorEnd = g4eTrajState.GetError();
0387
0388   CurvilinearTrajectoryError curvError(
0389       TrackPropagation::g4ErrorTrajErrToAlgebraicSymMatrix55(g4errorEnd, ftsStart.charge()));
0390
0391   if (mode == G4ErrorMode_PropBackwards) {
0392     GlobalTrajectoryParameters endParm(
0393         posEndGV, momEndGV, ftsStart.parameters().charge(), &ftsStart.parameters().magneticField());
0394
0395     // flip the momentum direction because it has been flipped before running
0396     // G4's backwards prop
0397     momEndGV = -momEndGV;
0398   }
0399
0400   LogDebug("Geant4e") << "G4e -  Error matrix after propagation: " << std::endl << g4errorEnd;
0401
0402   LogDebug("Geant4e") << "CMS -  Error matrix after propagation: " << std::endl << curvError.matrix();
0403
0404   GlobalTrajectoryParameters tParsDest(posEndGV, momEndGV, ftsStart.charge(), theField);
0405
0406   SurfaceSideDefinition::SurfaceSide side;
0407
0408   side = propagationDirection() == alongMomentum ? SurfaceSideDefinition::afterSurface
0409                                                  : SurfaceSideDefinition::beforeSurface;
0410
0411   return TsosPP(TrajectoryStateOnSurface(tParsDest, curvError, pDest, side), finalPathLength);
0412 }
0413
0414 //
0415 ////////////////////////////////////////////////////////////////////////////
0416 //
0417
0418 /** The methods propagateWithPath() are identical to the corresponding
0419  *  methods propagate() in what concerns the resulting
0420  *  TrajectoryStateOnSurface, but they provide in addition the
0421  *  exact path length along the trajectory.
0422  */
0423
0424 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
0425                                                                                  const Plane &pDest) const {
0426   // Finally build the pair<...> that needs to be returned where the second
0427   // parameter is the exact path length. Currently calculated with a stepping
0428   // action that adds up the length of every step
0429   return propagateGeneric(ftsStart, pDest);
0430 }
0431
0432 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
0433                                                                                  const Cylinder &cDest) const {
0434   // Finally build the pair<...> that needs to be returned where the second
0435   // parameter is the exact path length.
0436   return propagateGeneric(ftsStart, cDest);
0437 }
0438
0439 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(
0440     const TrajectoryStateOnSurface &tsosStart, const Plane &pDest) const {
0441   // Finally build the pair<...> that needs to be returned where the second
0442   // parameter is the exact path length.
0443   const FreeTrajectoryState ftsStart = *tsosStart.freeState();
0444   return propagateGeneric(ftsStart, pDest);
0445 }
0446
0447 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(
0448     const TrajectoryStateOnSurface &tsosStart, const Cylinder &cDest) const {
0449   const FreeTrajectoryState ftsStart = *tsosStart.freeState();
0450   // Finally build the pair<...> that needs to be returned where the second
0451   // parameter is the exact path length.
0452   return propagateGeneric(ftsStart, cDest);
0453 }
0454
0455 void Geant4ePropagator::debugReportPlaneSetup(GlobalPoint const &posPlane,
0456                                               HepGeom::Point3D<double> const &surfPos,
0457                                               GlobalVector const &normalPlane,
0458                                               HepGeom::Normal3D<double> const &surfNorm,
0459                                               const Plane &pDest) const {
0460   LogDebug("Geant4e") << "G4e -  Destination CMS plane position:" << posPlane << "cm\n"
0461                       << "G4e -                  (Ro, eta, phi): (" << posPlane.perp() << " cm, " << posPlane.eta()
0462                       << ", " << posPlane.phi().degrees() << " deg)\n"
0463                       << "G4e -  Destination G4  plane position: " << surfPos << " mm, Ro = " << surfPos.perp()
0464                       << " mm";
0465   LogDebug("Geant4e") << "G4e -  Destination CMS plane normal  : " << normalPlane << "\n"
0466                       << "G4e -  Destination G4  plane normal  : " << normalPlane;
0467   LogDebug("Geant4e") << "G4e -  Distance from plane position to plane: " << pDest.localZ(posPlane) << " cm";
0468 }
0469
0470 template <class SurfaceType>
0471 void Geant4ePropagator::debugReportTrackState(std::string const &currentContext,
0472                                               GlobalPoint const &cmsInitPos,
0473                                               CLHEP::Hep3Vector const &g4InitPos,
0474                                               GlobalVector const &cmsInitMom,
0475                                               CLHEP::Hep3Vector const &g4InitMom,
0476                                               const SurfaceType &pDest) const {
0477   LogDebug("Geant4e") << "G4e - Current Context: " << currentContext;
0478   LogDebug("Geant4e") << "G4e -  CMS point position:" << cmsInitPos << "cm\n"
0479                       << "G4e -              (Ro, eta, phi): (" << cmsInitPos.perp() << " cm, " << cmsInitPos.eta()
0480                       << ", " << cmsInitPos.phi().degrees() << " deg)\n"
0481                       << "G4e -   G4  point position: " << g4InitPos << " mm, Ro = " << g4InitPos.perp() << " mm";
0482   LogDebug("Geant4e") << "G4e -   CMS momentum      :" << cmsInitMom << "GeV\n"
0483                       << " pt: " << cmsInitMom.perp() << "G4e -  G4  momentum      : " << g4InitMom << " MeV";
0484 }