Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-01-21 00:19:57

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