File indexing completed on 2024-05-10 02:21:31
0001 #include <sstream>
0002
0003
0004 #include "TrackPropagation/Geant4e/interface/ConvertFromToCLHEP.h"
0005 #include "TrackPropagation/Geant4e/interface/Geant4ePropagator.h"
0006
0007
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
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
0038 #include <CLHEP/Units/SystemOfUnits.h>
0039
0040
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
0055
0056 ensureGeant4eIsInitilized(true);
0057 }
0058
0059
0060
0061 Geant4ePropagator::~Geant4ePropagator() {
0062 LogDebug("Geant4e") << "Geant4ePropagator::~Geant4ePropagator()" << std::endl;
0063
0064
0065
0066
0067 }
0068
0069
0070
0071
0072
0073
0074
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
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
0099 GlobalPoint posPlane = pDest.toGlobal(LocalPoint(0, 0, 0));
0100 GlobalVector normalPlane = pDest.toGlobal(LocalVector(0, 0, 1.));
0101 normalPlane = normalPlane.unit();
0102
0103
0104
0105
0106
0107 HepGeom::Point3D<double> surfPos = TrackPropagation::globalPointToHepPoint3D(posPlane);
0108 HepGeom::Normal3D<double> surfNorm = TrackPropagation::globalVectorToHepNormal3D(normalPlane);
0109
0110
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
0118
0119
0120 G4float radCyl = pDest.radius() * CLHEP::cm;
0121
0122 G4ThreeVector posCyl = TrackPropagation::globalPointToHep3Vector(pDest.position());
0123
0124 G4RotationMatrix rotCyl = TrackPropagation::tkRotationFToHepRotation(pDest.rotation());
0125
0126
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
0190
0191 LocalPoint lpos = pDest.toLocal(cmsInitPos);
0192 Surface::Side theSide = pDest.side(lpos, 0);
0193 if (theSide == SurfaceOrientation::positiveSide) {
0194 mode = G4ErrorMode_PropBackwards;
0195 LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\' indirect "
0196 "via the Any direction";
0197 } else {
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
0234
0235
0236
0237 ErrorTargetPair g4eTarget_center = transformToG4SurfaceTarget(pDest, false);
0238
0239
0240
0241
0242 GlobalPoint cmsInitPos = ftsStart.position();
0243 GlobalVector cmsInitMom = ftsStart.momentum();
0244 bool flipped = false;
0245 if (propagationDirection() == oppositeToMomentum) {
0246
0247
0248 cmsInitMom = -cmsInitMom;
0249 flipped = true;
0250 }
0251
0252
0253 G4ErrorMode mode = G4ErrorMode_PropForwards;
0254 if (!configurePropagation(mode, pDest, cmsInitPos, cmsInitMom))
0255 return TsosPP(TrajectoryStateOnSurface(), 0.0f);
0256
0257
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
0267
0268
0269
0270
0271
0272
0273
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
0289
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
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
0314 theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointAndSetup(
0315 g4InitPos, &g4InitMom, false, false);
0316
0317 bool continuePropagation = true;
0318 while (continuePropagation) {
0319 iterations++;
0320 LogDebug("Geant4e") << std::endl << "step count " << iterations << " step length " << finalPathLength;
0321
0322
0323 theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointWithinVolume(g4eTrajState.GetPosition());
0324
0325 const int ierr = theG4eManager->PropagateOneStep(&g4eTrajState, mode);
0326
0327 if (ierr != 0) {
0328
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
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
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
0359 if (propagationDirection() == oppositeToMomentum)
0360 finalPathLength = -finalPathLength;
0361
0362
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
0372
0373
0374
0375 const HepGeom::Vector3D<double> momEnd = g4eTrajState.GetMomentum();
0376
0377
0378
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
0385
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
0396
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
0419
0420
0421
0422
0423
0424 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
0425 const Plane &pDest) const {
0426
0427
0428
0429 return propagateGeneric(ftsStart, pDest);
0430 }
0431
0432 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
0433 const Cylinder &cDest) const {
0434
0435
0436 return propagateGeneric(ftsStart, cDest);
0437 }
0438
0439 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(
0440 const TrajectoryStateOnSurface &tsosStart, const Plane &pDest) const {
0441
0442
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
0451
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 ¤tContext,
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 }