File indexing completed on 2023-01-21 00:19:57
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
0036
0037 #include "CLHEP/Units/GlobalSystemOfUnits.h"
0038
0039
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
0054
0055 ensureGeant4eIsInitilized(true);
0056 }
0057
0058
0059
0060 Geant4ePropagator::~Geant4ePropagator() {
0061 LogDebug("Geant4e") << "Geant4ePropagator::~Geant4ePropagator()" << std::endl;
0062
0063
0064
0065
0066 }
0067
0068
0069
0070
0071
0072
0073
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
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
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
0101 GlobalPoint posPlane = pDest.toGlobal(LocalPoint(0, 0, 0));
0102 GlobalVector normalPlane = pDest.toGlobal(LocalVector(0, 0, 1.));
0103 normalPlane = normalPlane.unit();
0104
0105
0106
0107
0108
0109 HepGeom::Point3D<double> surfPos = TrackPropagation::globalPointToHepPoint3D(posPlane);
0110 HepGeom::Normal3D<double> surfNorm = TrackPropagation::globalVectorToHepNormal3D(normalPlane);
0111
0112
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
0120
0121
0122 G4float radCyl = pDest.radius() * cm;
0123
0124 G4ThreeVector posCyl = TrackPropagation::globalPointToHep3Vector(pDest.position());
0125
0126 G4RotationMatrix rotCyl = TrackPropagation::tkRotationFToHepRotation(pDest.rotation());
0127
0128
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
0192
0193 LocalPoint lpos = pDest.toLocal(cmsInitPos);
0194 Surface::Side theSide = pDest.side(lpos, 0);
0195 if (theSide == SurfaceOrientation::positiveSide) {
0196 mode = G4ErrorMode_PropBackwards;
0197 LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\' indirect "
0198 "via the Any direction";
0199 } else {
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
0236
0237
0238
0239 ErrorTargetPair g4eTarget_center = transformToG4SurfaceTarget(pDest, false);
0240
0241
0242
0243
0244 GlobalPoint cmsInitPos = ftsStart.position();
0245 GlobalVector cmsInitMom = ftsStart.momentum();
0246 bool flipped = false;
0247 if (propagationDirection() == oppositeToMomentum) {
0248
0249
0250 cmsInitMom = -cmsInitMom;
0251 flipped = true;
0252 }
0253
0254
0255 G4ErrorMode mode = G4ErrorMode_PropForwards;
0256 if (!configurePropagation(mode, pDest, cmsInitPos, cmsInitMom))
0257 return TsosPP(TrajectoryStateOnSurface(), 0.0f);
0258
0259
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
0269
0270
0271
0272
0273
0274
0275
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
0291
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
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
0316 theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointAndSetup(
0317 g4InitPos, &g4InitMom, false, false);
0318
0319 bool continuePropagation = true;
0320 while (continuePropagation) {
0321 iterations++;
0322 LogDebug("Geant4e") << std::endl << "step count " << iterations << " step length " << finalPathLength;
0323
0324
0325 theG4eManager->GetErrorPropagationNavigator()->LocateGlobalPointWithinVolume(g4eTrajState.GetPosition());
0326
0327 const int ierr = theG4eManager->PropagateOneStep(&g4eTrajState, mode);
0328
0329 if (ierr != 0) {
0330
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
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
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
0361 if (propagationDirection() == oppositeToMomentum)
0362 finalPathLength = -finalPathLength;
0363
0364
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
0374
0375
0376
0377 const HepGeom::Vector3D<double> momEnd = g4eTrajState.GetMomentum();
0378
0379
0380
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
0387
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
0398
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
0421
0422
0423
0424
0425
0426 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
0427 const Plane &pDest) const {
0428
0429
0430
0431 return propagateGeneric(ftsStart, pDest);
0432 }
0433
0434 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(const FreeTrajectoryState &ftsStart,
0435 const Cylinder &cDest) const {
0436
0437
0438 return propagateGeneric(ftsStart, cDest);
0439 }
0440
0441 std::pair<TrajectoryStateOnSurface, double> Geant4ePropagator::propagateWithPath(
0442 const TrajectoryStateOnSurface &tsosStart, const Plane &pDest) const {
0443
0444
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
0453
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 ¤tContext,
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 }