CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Geant4ePropagator.cc
Go to the documentation of this file.
1 #include <sstream>
2 
3 // Geant4e
6 
7 // CMSSW
12 
17 
18 // Geant4
19 #include "G4Box.hh"
20 #include "G4ErrorCylSurfaceTarget.hh"
21 #include "G4ErrorFreeTrajState.hh"
22 #include "G4ErrorPlaneSurfaceTarget.hh"
23 #include "G4ErrorPropagatorData.hh"
24 #include "G4ErrorRunManagerHelper.hh"
25 #include "G4EventManager.hh"
26 #include "G4Field.hh"
27 #include "G4FieldManager.hh"
28 #include "G4GeometryTolerance.hh"
29 #include "G4SteppingControl.hh"
30 #include "G4TransportationManager.hh"
31 #include "G4Tubs.hh"
32 #include "G4UImanager.hh"
33 #include "G4ErrorPropagationNavigator.hh"
34 #include "G4RunManagerKernel.hh"
35 #include "G4StateManager.hh"
36 
37 // CLHEP
38 #include "CLHEP/Units/GlobalSystemOfUnits.h"
39 
45  double plimit)
46  : Propagator(dir),
47  theField(field),
48  theParticleName(particleName),
49  theG4eManager(G4ErrorPropagatorManager::GetErrorPropagatorManager()),
50  theG4eData(G4ErrorPropagatorData::GetErrorPropagatorData()),
51  plimit_(plimit) {
52  LogDebug("Geant4e") << "Geant4e Propagator initialized";
53 
54  // has to be called here, doing it later will not load the G4 physics list
55  // properly when using the G4 ES Producer. Reason: unclear
57 }
58 
62  LogDebug("Geant4e") << "Geant4ePropagator::~Geant4ePropagator()" << std::endl;
63 
64  // don't close the g4 Geometry here, because the propagator might have been
65  // cloned
66  // but there is only one, globally-shared Geometry
67 }
68 
69 //
71 //
72 
78  LogDebug("Geant4ePropagator") << "G4 propagator starts isInitialized, theField: " << theField;
79 
80  auto man = G4RunManagerKernel::GetRunManagerKernel();
81  if (G4StateManager::GetStateManager()->GetCurrentState() == G4State_PreInit) {
82  man->SetVerboseLevel(0);
83  theG4eManager->InitGeant4e();
84 
85  // define 10 mm step limit for propagator
86  G4UImanager::GetUIpointer()->ApplyCommand("/geant4e/limits/stepLength 10.0 mm");
87  }
88  const G4Field *field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
89  if (field == nullptr) {
90  edm::LogError("Geant4e") << "No G4 magnetic field defined";
91  }
92  LogDebug("Geant4ePropagator") << "G4 propagator initialized; field: " << field;
93 }
94 
95 template <>
97  bool moveTargetToEndOfSurface) const {
98  //* Get position and normal (orientation) of the destination plane
99  GlobalPoint posPlane = pDest.toGlobal(LocalPoint(0, 0, 0));
100  GlobalVector normalPlane = pDest.toGlobal(LocalVector(0, 0, 1.));
101  normalPlane = normalPlane.unit();
102 
103  //* Transform this into HepGeom::Point3D<double> and
104  // HepGeom::Normal3D<double> that define a plane for
105  // Geant4e.
106  // CMS uses cm and GeV while Geant4 uses mm and MeV
107  HepGeom::Point3D<double> surfPos = TrackPropagation::globalPointToHepPoint3D(posPlane);
108  HepGeom::Normal3D<double> surfNorm = TrackPropagation::globalVectorToHepNormal3D(normalPlane);
109 
110  //* Set the target surface
111  return ErrorTargetPair(false, std::make_shared<G4ErrorPlaneSurfaceTarget>(surfNorm, surfPos));
112 }
113 
114 template <>
116  bool moveTargetToEndOfSurface) const {
117  // Get Cylinder parameters.
118  // CMS uses cm and GeV while Geant4 uses mm and MeV.
119  // - Radius
120  G4float radCyl = pDest.radius() * cm;
121  // - Position: PositionType & GlobalPoint are Basic3DPoint<float,GlobalTag>
122  G4ThreeVector posCyl = TrackPropagation::globalPointToHep3Vector(pDest.position());
123  // - Rotation: Type in CMSSW is RotationType == TkRotation<T>, T=float
124  G4RotationMatrix rotCyl = TrackPropagation::tkRotationFToHepRotation(pDest.rotation());
125 
126  // DEBUG
128  LogDebug("Geant4e") << "G4e - TkRotation" << rotation;
129  LogDebug("Geant4e") << "G4e - G4Rotation" << rotCyl << "mm";
130 
131  return ErrorTargetPair(!moveTargetToEndOfSurface, std::make_shared<G4ErrorCylSurfaceTarget>(radCyl, posCyl, rotCyl));
132 }
133 
134 template <>
136  return "Cylinder";
137 }
138 
139 template <>
141  return "Plane";
142 }
143 
146 
147  if (charge > 0) {
148  particleName += "+";
149  }
150  if (charge < 0) {
151  particleName += "-";
152  }
153 
154  LogDebug("Geant4e") << "G4e - Particle name: " << particleName;
155 
156  return particleName;
157 }
158 
159 template <>
161  Plane const &pDest,
162  GlobalPoint const &cmsInitPos,
163  GlobalVector const &cmsInitMom) const {
164  if (cmsInitMom.mag() < plimit_)
165  return false;
166  if (pDest.localZ(cmsInitPos) * pDest.localZ(cmsInitMom) < 0) {
167  mode = G4ErrorMode_PropForwards;
168  LogDebug("Geant4e") << "G4e - Propagator mode is \'forwards\' indirect "
169  "via the Any direction"
170  << std::endl;
171  } else {
172  mode = G4ErrorMode_PropBackwards;
173  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\' indirect "
174  "via the Any direction"
175  << std::endl;
176  }
177 
178  return true;
179 }
180 
181 template <>
183  Cylinder const &pDest,
184  GlobalPoint const &cmsInitPos,
185  GlobalVector const &cmsInitMom) const {
186  if (cmsInitMom.mag() < plimit_)
187  return false;
188  //------------------------------------
189  // For cylinder assume outside is backwards, inside is along
190  // General use for particles from collisions
191  LocalPoint lpos = pDest.toLocal(cmsInitPos);
192  Surface::Side theSide = pDest.side(lpos, 0);
193  if (theSide == SurfaceOrientation::positiveSide) { // outside cylinder
194  mode = G4ErrorMode_PropBackwards;
195  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\' indirect "
196  "via the Any direction";
197  } else { // inside cylinder
198  mode = G4ErrorMode_PropForwards;
199  LogDebug("Geant4e") << "G4e - Propagator mode is \'forwards\' indirect "
200  "via the Any direction";
201  }
202 
203  return true;
204 }
205 
206 template <class SurfaceType>
208  SurfaceType const &pDest,
209  GlobalPoint const &cmsInitPos,
210  GlobalVector const &cmsInitMom) const {
211  if (cmsInitMom.mag() < plimit_)
212  return false;
214  mode = G4ErrorMode_PropBackwards;
215  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\' " << std::endl;
216  } else if (propagationDirection() == alongMomentum) {
217  mode = G4ErrorMode_PropForwards;
218  LogDebug("Geant4e") << "G4e - Propagator mode is \'forwards\'" << std::endl;
219  } else if (propagationDirection() == anyDirection) {
220  if (configureAnyPropagation(mode, pDest, cmsInitPos, cmsInitMom) == false)
221  return false;
222  } else {
223  edm::LogError("Geant4e") << "G4e - Unsupported propagation mode";
224  return false;
225  }
226  return true;
227 }
228 
229 template <class SurfaceType>
230 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateGeneric(const FreeTrajectoryState &ftsStart,
231  const SurfaceType &pDest) const {
233  // Construct the target surface
234  //
235  //* Set the target surface
236 
237  ErrorTargetPair g4eTarget_center = transformToG4SurfaceTarget(pDest, false);
238 
239  // * Get the starting point and direction and convert them to
240  // CLHEP::Hep3Vector
241  // for G4. CMS uses cm and GeV while Geant4 uses mm and MeV
242  GlobalPoint cmsInitPos = ftsStart.position();
243  GlobalVector cmsInitMom = ftsStart.momentum();
244  bool flipped = false;
246  // flip the momentum vector as Geant4 will not do this
247  // on it's own in a backward propagation
248  cmsInitMom = -cmsInitMom;
249  flipped = true;
250  }
251 
252  // Set the mode of propagation according to the propagation direction
253  G4ErrorMode mode = G4ErrorMode_PropForwards;
254  if (!configurePropagation(mode, pDest, cmsInitPos, cmsInitMom))
255  return TsosPP(TrajectoryStateOnSurface(), 0.0f);
256 
257  // re-check propagation direction chosen in case of AnyDirection
258  if (mode == G4ErrorMode_PropBackwards && !flipped)
259  cmsInitMom = -cmsInitMom;
260 
261  CLHEP::Hep3Vector g4InitPos = TrackPropagation::globalPointToHep3Vector(cmsInitPos);
262  CLHEP::Hep3Vector g4InitMom = TrackPropagation::globalVectorToHep3Vector(cmsInitMom * GeV);
263 
264  debugReportTrackState("intitial", cmsInitPos, g4InitPos, cmsInitMom, g4InitMom, pDest);
265 
266  // Set the mode of propagation according to the propagation direction
267  // G4ErrorMode mode = G4ErrorMode_PropForwards;
268 
269  // if (!configurePropagation(mode, pDest, cmsInitPos, cmsInitMom))
270  // return TsosPP(TrajectoryStateOnSurface(), 0.0f);
271 
273  // Set the error and trajectories, and finally propagate
274  //
275  G4ErrorTrajErr g4error(5, 1);
276  if (ftsStart.hasError()) {
278  initErr = ftsStart.curvilinearError();
279  g4error = TrackPropagation::algebraicSymMatrix55ToG4ErrorTrajErr(initErr, ftsStart.charge());
280  LogDebug("Geant4e") << "CMS - Error matrix: " << std::endl << initErr.matrix();
281  } else {
282  LogDebug("Geant4e") << "No error matrix available" << std::endl;
283  return TsosPP(TrajectoryStateOnSurface(), 0.0f);
284  }
285 
286  LogDebug("Geant4e") << "G4e - Error matrix: " << std::endl << g4error;
287 
288  // in CMSSW, the state errors are deflated when performing the backward
289  // propagation
290  if (mode == G4ErrorMode_PropForwards) {
291  G4ErrorPropagatorData::GetErrorPropagatorData()->SetStage(G4ErrorStage_Inflation);
292  } else if (mode == G4ErrorMode_PropBackwards) {
293  G4ErrorPropagatorData::GetErrorPropagatorData()->SetStage(G4ErrorStage_Deflation);
294  }
295 
296  G4ErrorFreeTrajState g4eTrajState(generateParticleName(ftsStart.charge()), g4InitPos, g4InitMom, g4error);
297  LogDebug("Geant4e") << "G4e - Traj. State: " << (g4eTrajState);
298 
300  // Propagate
301  int iterations = 0;
302  double finalPathLength = 0;
303 
304  HepGeom::Point3D<double> finalRecoPos;
305 
306  G4ErrorPropagatorData::GetErrorPropagatorData()->SetMode(mode);
307 
308  theG4eData->SetTarget(g4eTarget_center.second.get());
309  LogDebug("Geant4e") << "Running Propagation to the RECO surface" << std::endl;
310 
311  theG4eManager->InitTrackPropagation();
312 
313  // re-initialize navigator to avoid mismatches and/or segfaults
314  theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointAndSetup(
315  g4InitPos, &g4InitMom, /*pRelativeSearch = */ false, /*ignoreDirection = */ false);
316 
317  bool continuePropagation = true;
318  while (continuePropagation) {
319  iterations++;
320  LogDebug("Geant4e") << std::endl << "step count " << iterations << " step length " << finalPathLength;
321 
322  // re-initialize navigator to avoid mismatches and/or segfaults
323  theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointWithinVolume(g4eTrajState.GetPosition());
324 
325  const int ierr = theG4eManager->PropagateOneStep(&g4eTrajState, mode);
326 
327  if (ierr != 0) {
328  // propagation failed, return invalid track state
329  return TsosPP(TrajectoryStateOnSurface(), 0.0f);
330  }
331 
332  const float thisPathLength = TrackPropagation::g4doubleToCmsDouble(g4eTrajState.GetG4Track()->GetStepLength());
333 
334  LogDebug("Geant4e") << "step Length was " << thisPathLength << " cm, current global position: "
335  << TrackPropagation::hepPoint3DToGlobalPoint(g4eTrajState.GetPosition()) << std::endl;
336 
337  finalPathLength += thisPathLength;
338 
339  // if (std::fabs(finalPathLength) > 10000.0f)
340  if (std::fabs(finalPathLength) > 200.0f) {
341  LogDebug("Geant4e") << "ERROR: Quitting propagation: path length mega large" << std::endl;
342  theG4eManager->GetPropagator()->InvokePostUserTrackingAction(g4eTrajState.GetG4Track());
343  continuePropagation = false;
344  LogDebug("Geant4e") << "WARNING: Quitting propagation: max path length "
345  "exceeded, returning invalid state"
346  << std::endl;
347 
348  // reached maximum path length, bail out
349  return TsosPP(TrajectoryStateOnSurface(), 0.0f);
350  }
351 
352  if (theG4eManager->GetPropagator()->CheckIfLastStep(g4eTrajState.GetG4Track())) {
353  theG4eManager->GetPropagator()->InvokePostUserTrackingAction(g4eTrajState.GetG4Track());
354  continuePropagation = false;
355  }
356  }
357 
358  // CMSSW Tracking convention, backward propagations have negative path length
360  finalPathLength = -finalPathLength;
361 
362  // store the correct location for the hit on the RECO surface
363  LogDebug("Geant4e") << "Position on the RECO surface" << g4eTrajState.GetPosition() << std::endl;
364  finalRecoPos = g4eTrajState.GetPosition();
365 
366  theG4eManager->EventTermination();
367 
368  LogDebug("Geant4e") << "Final position of the Track :" << g4eTrajState.GetPosition() << std::endl;
369 
371  // Retrieve the state in the end from Geant4e, convert them to CMS vectors
372  // and points, and build global trajectory parameters.
373  // CMS uses cm and GeV while Geant4 uses mm and MeV
374  //
375  const HepGeom::Vector3D<double> momEnd = g4eTrajState.GetMomentum();
376 
377  // use the hit on the the RECO plane as the final position to be d'accor with
378  // the RecHit measurements
379  const GlobalPoint posEndGV = TrackPropagation::hepPoint3DToGlobalPoint(finalRecoPos);
381 
382  debugReportTrackState("final", posEndGV, finalRecoPos, momEndGV, momEnd, pDest);
383 
384  // Get the error covariance matrix from Geant4e. It comes in curvilinear
385  // coordinates so use the appropiate CMS class
386  G4ErrorTrajErr g4errorEnd = g4eTrajState.GetError();
387 
388  CurvilinearTrajectoryError curvError(
390 
391  if (mode == G4ErrorMode_PropBackwards) {
393  posEndGV, momEndGV, ftsStart.parameters().charge(), &ftsStart.parameters().magneticField());
394 
395  // flip the momentum direction because it has been flipped before running
396  // G4's backwards prop
397  momEndGV = -momEndGV;
398  }
399 
400  LogDebug("Geant4e") << "G4e - Error matrix after propagation: " << std::endl << g4errorEnd;
401 
402  LogDebug("Geant4e") << "CMS - Error matrix after propagation: " << std::endl << curvError.matrix();
403 
404  GlobalTrajectoryParameters tParsDest(posEndGV, momEndGV, ftsStart.charge(), theField);
405 
407 
410 
411  return TsosPP(TrajectoryStateOnSurface(tParsDest, curvError, pDest, side), finalPathLength);
412 }
413 
414 //
416 //
417 
424 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
425  const Plane &pDest) const {
426  // Finally build the pair<...> that needs to be returned where the second
427  // parameter is the exact path length. Currently calculated with a stepping
428  // action that adds up the length of every step
429  return propagateGeneric(ftsStart, pDest);
430 }
431 
432 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
433  const Cylinder &cDest) const {
434  // Finally build the pair<...> that needs to be returned where the second
435  // parameter is the exact path length.
436  return propagateGeneric(ftsStart, cDest);
437 }
438 
439 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(
440  const TrajectoryStateOnSurface &tsosStart, const Plane &pDest) const {
441  // Finally build the pair<...> that needs to be returned where the second
442  // parameter is the exact path length.
443  const FreeTrajectoryState ftsStart = *tsosStart.freeState();
444  return propagateGeneric(ftsStart, pDest);
445 }
446 
447 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(
448  const TrajectoryStateOnSurface &tsosStart, const Cylinder &cDest) const {
449  const FreeTrajectoryState ftsStart = *tsosStart.freeState();
450  // Finally build the pair<...> that needs to be returned where the second
451  // parameter is the exact path length.
452  return propagateGeneric(ftsStart, cDest);
453 }
454 
456  HepGeom::Point3D<double> const &surfPos,
457  GlobalVector const &normalPlane,
458  HepGeom::Normal3D<double> const &surfNorm,
459  const Plane &pDest) const {
460  LogDebug("Geant4e") << "G4e - Destination CMS plane position:" << posPlane << "cm\n"
461  << "G4e - (Ro, eta, phi): (" << posPlane.perp() << " cm, " << posPlane.eta()
462  << ", " << posPlane.phi().degrees() << " deg)\n"
463  << "G4e - Destination G4 plane position: " << surfPos << " mm, Ro = " << surfPos.perp()
464  << " mm";
465  LogDebug("Geant4e") << "G4e - Destination CMS plane normal : " << normalPlane << "\n"
466  << "G4e - Destination G4 plane normal : " << normalPlane;
467  LogDebug("Geant4e") << "G4e - Distance from plane position to plane: " << pDest.localZ(posPlane) << " cm";
468 }
469 
470 template <class SurfaceType>
472  GlobalPoint const &cmsInitPos,
473  CLHEP::Hep3Vector const &g4InitPos,
474  GlobalVector const &cmsInitMom,
475  CLHEP::Hep3Vector const &g4InitMom,
476  const SurfaceType &pDest) const {
477  LogDebug("Geant4e") << "G4e - Current Context: " << currentContext;
478  LogDebug("Geant4e") << "G4e - CMS point position:" << cmsInitPos << "cm\n"
479  << "G4e - (Ro, eta, phi): (" << cmsInitPos.perp() << " cm, " << cmsInitPos.eta()
480  << ", " << cmsInitPos.phi().degrees() << " deg)\n"
481  << "G4e - G4 point position: " << g4InitPos << " mm, Ro = " << g4InitPos.perp() << " mm";
482  LogDebug("Geant4e") << "G4e - CMS momentum :" << cmsInitMom << "GeV\n"
483  << " pt: " << cmsInitMom.perp() << "G4e - G4 momentum : " << g4InitMom << " MeV";
484 }
std::pair< TrajectoryStateOnSurface, double > TsosPP
Local3DVector LocalVector
Definition: LocalVector.h:12
T perp() const
Definition: PV3DBase.h:69
Point3DBase< Scalar, LocalTag > LocalPoint
Definition: Definitions.h:30
const CurvilinearTrajectoryError & curvilinearError() const
void debugReportTrackState(std::string const &currentContext, GlobalPoint const &cmsInitPos, CLHEP::Hep3Vector const &g4InitPos, GlobalVector const &cmsInitMom, CLHEP::Hep3Vector const &g4InitMom, const SurfaceType &pDest) const
HepGeom::Normal3D< double > globalVectorToHepNormal3D(const GlobalVector &p)
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
double g4doubleToCmsDouble(const G4double &d)
T eta() const
Definition: PV3DBase.h:73
std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &, const Plane &) const override
std::string theParticleName
virtual PropagationDirection propagationDirection() const final
Definition: Propagator.h:139
CLHEP::Hep3Vector globalVectorToHep3Vector(const GlobalVector &p)
PropagationDirection
Log< level::Error, false > LogError
const MagneticField * theField
LocalPoint toLocal(const GlobalPoint &gp) const
const GlobalTrajectoryParameters & parameters() const
Definition: Plane.h:16
G4ErrorPropagatorManager * theG4eManager
void ensureGeant4eIsInitilized(bool forceInit) const
GlobalPoint position() const
CLHEP::HepRotation tkRotationFToHepRotation(const TkRotation< float > &tkr)
AlgebraicSymMatrix55 g4ErrorTrajErrToAlgebraicSymMatrix55(const G4ErrorTrajErr &e, const int q)
G4ErrorPropagatorData * theG4eData
TrackCharge charge() const
GlobalPoint hepPoint3DToGlobalPoint(const HepGeom::Point3D< double > &r)
std::pair< bool, std::shared_ptr< G4ErrorTarget > > ErrorTargetPair
T mag() const
Definition: PV3DBase.h:64
GlobalVector momentum() const
double f[11][100]
void debugReportPlaneSetup(GlobalPoint const &posPlane, HepGeom::Point3D< double > const &surfPos, GlobalVector const &normalPlane, HepGeom::Normal3D< double > const &surfNorm, const Plane &pDest) const
bool configureAnyPropagation(G4ErrorMode &mode, SurfaceType const &pDest, GlobalPoint const &cmsInitPos, GlobalVector const &cmsInitMom) const
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:79
float localZ(const GlobalPoint &gp) const
Definition: Plane.h:45
HepGeom::Point3D< double > globalPointToHepPoint3D(const GlobalPoint &r)
const PositionType & position() const
~Geant4ePropagator() override
const AlgebraicSymMatrix55 & matrix() const
T1 degrees() const
Definition: Phi.h:107
const MagneticField & magneticField() const
Side side(const LocalPoint &p, Scalar toler) const override
Definition: Cylinder.cc:9
GlobalVector hep3VectorToGlobalVector(const CLHEP::Hep3Vector &p)
std::string getSurfaceType(SurfaceType const &surface) const
std::string generateParticleName(int charge) const
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:64
FreeTrajectoryState const * freeState(bool withErrors=true) const
bool configurePropagation(G4ErrorMode &mode, SurfaceType const &pDest, GlobalPoint const &cmsInitPos, GlobalVector const &cmsInitMom) const
G4ErrorTrajErr algebraicSymMatrix55ToG4ErrorTrajErr(const AlgebraicSymMatrix55 &e, const int q)
const RotationType & rotation() const
Vector3DBase unit() const
Definition: Vector3DBase.h:54
ErrorTargetPair transformToG4SurfaceTarget(const SurfaceType &pDest, bool moveTargetToEndOfSurface) const
Geant4ePropagator(const MagneticField *field=nullptr, std::string particleName="mu", PropagationDirection dir=alongMomentum, double plimit=1.0)
CLHEP::Hep3Vector globalPointToHep3Vector(const GlobalPoint &r)
std::pair< TrajectoryStateOnSurface, double > propagateGeneric(const FreeTrajectoryState &ftsStart, const SurfaceType &pDest) const
#define LogDebug(id)