17 #include "G4ErrorFreeTrajState.hh"
18 #include "G4ErrorPlaneSurfaceTarget.hh"
19 #include "G4ErrorCylSurfaceTarget.hh"
20 #include "G4ErrorPropagatorData.hh"
21 #include "G4EventManager.hh"
22 #include "G4SteppingControl.hh"
25 #include "CLHEP/Units/GlobalSystemOfUnits.h"
31 const char* particleName,
35 theParticleName(particleName),
36 theG4eManager(G4ErrorPropagatorManager::GetErrorPropagatorManager()),
37 theSteppingAction(0) {
39 G4ErrorPropagatorData::SetVerbose(0);
44 Geant4ePropagator::~Geant4ePropagator() {
57 const Plane& pDest)
const {
59 if(theG4eManager->PrintG4ErrorState() ==
"G4ErrorState_PreInit")
60 theG4eManager->InitGeant4e();
62 if (!theSteppingAction) {
63 theSteppingAction =
new Geant4eSteppingAction;
64 theG4eManager->SetUserAction(theSteppingAction);
74 normalPlane = normalPlane.
unit();
79 HepGeom::Point3D<double> surfPos =
81 HepGeom::Normal3D<double> surfNorm =
85 LogDebug(
"Geant4e") <<
"G4e - Destination CMS plane position:" << posPlane <<
"cm\n"
86 <<
"G4e - (Ro, eta, phi): ("
87 << posPlane.
perp() <<
" cm, "
88 << posPlane.
eta() <<
", "
89 << posPlane.
phi().degrees() <<
" deg)\n"
90 <<
"G4e - Destination G4 plane position: " << surfPos
91 <<
" mm, Ro = " << surfPos.perp() <<
" mm";
92 LogDebug(
"Geant4e") <<
"G4e - Destination CMS plane normal : "
93 << normalPlane <<
"\n"
94 <<
"G4e - Destination G4 plane normal : "
96 LogDebug(
"Geant4e") <<
"G4e - Distance from plane position to plane: "
97 << pDest.
localZ(posPlane) <<
" cm";
101 G4ErrorSurfaceTarget* g4eTarget =
new G4ErrorPlaneSurfaceTarget(surfNorm,
117 CLHEP::Hep3Vector g4InitPos =
119 CLHEP::Hep3Vector g4InitMom =
123 LogDebug(
"Geant4e") <<
"G4e - Initial CMS point position:" << cmsInitPos
125 <<
"G4e - (Ro, eta, phi): ("
126 << cmsInitPos.
perp() <<
" cm, "
127 << cmsInitPos.
eta() <<
", "
128 << cmsInitPos.
phi().degrees() <<
" deg)\n"
129 <<
"G4e - Initial G4 point position: " << g4InitPos
130 <<
" mm, Ro = " << g4InitPos.perp() <<
" mm";
131 LogDebug(
"Geant4e") <<
"G4e - Initial CMS momentum :" << cmsInitMom
133 <<
"G4e - Initial G4 momentum : " << g4InitMom
135 LogDebug(
"Geant4e") <<
"G4e - Distance from initial point to plane: "
136 << pDest.
localZ(cmsInitPos) <<
" cm";
154 LogDebug(
"Geant4e") <<
"G4e - Particle name: " << particleName;
162 G4ErrorTrajErr g4error( 5, 1 );
167 LogDebug(
"Geant4e") <<
"G4e - Error matrix: " << g4error;
169 G4ErrorFreeTrajState* g4eTrajState =
170 new G4ErrorFreeTrajState(particleName, g4InitPos, g4InitMom, g4error);
171 LogDebug(
"Geant4e") <<
"G4e - Traj. State: " << (*g4eTrajState);
174 G4ErrorMode
mode = G4ErrorMode_PropForwards;
177 mode = G4ErrorMode_PropBackwards;
178 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'backwards\'";
180 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'forwards\'";
182 std::cout <<
"Determining actual direction";
183 if(pDest.
localZ(cmsInitPos)*pDest.
localZ(cmsInitMom) < 0) {
184 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'forwards\'";
185 std::cout <<
", got forwards" << std::endl;
187 mode = G4ErrorMode_PropBackwards;
188 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'backwards\'";
189 std::cout <<
", got backwards" << std::endl;
200 if(mode == G4ErrorMode_PropBackwards) {
204 g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
205 ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
206 g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
208 ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
210 LogDebug(
"Geant4e") <<
"G4e - Return error from propagation: " << ierr;
213 LogDebug(
"Geant4e") <<
"G4e - Error is not 0, returning invalid trajectory";
225 HepGeom::Point3D<double> posEnd = g4eTrajState->GetPosition();
226 HepGeom::Vector3D<double> momEnd = g4eTrajState->GetMomentum();
232 LogDebug(
"Geant4e") <<
"G4e - Final CMS point position:" << posEndGV
234 <<
"G4e - (Ro, eta, phi): ("
235 << posEndGV.
perp() <<
" cm, "
236 << posEndGV.
eta() <<
", "
237 << posEndGV.
phi().degrees() <<
" deg)\n"
238 <<
"G4e - Final G4 point position: " << posEnd
239 <<
" mm,\tRo =" << posEnd.perp() <<
" mm";
240 LogDebug(
"Geant4e") <<
"G4e - Final CMS momentum :" << momEndGV
242 <<
"G4e - Final G4 momentum : " << momEnd
244 LogDebug(
"Geant4e") <<
"G4e - Distance from final point to plane: "
245 << pDest.
localZ(posEndGV) <<
" cm";
253 G4ErrorTrajErr g4errorEnd = g4eTrajState->GetError();
256 LogDebug(
"Geant4e") <<
"G4e - Error matrix after propagation: " << g4errorEnd;
261 LogDebug(
"Geant4e") <<
"G4e - SurfaceSide is always atCenterOfSurface after propagation";
274 return propagate(ftsStart,plane);
285 if(theG4eManager->PrintG4ErrorState() ==
"G4ErrorState_PreInit")
286 theG4eManager->InitGeant4e();
287 if (!theSteppingAction) {
288 theSteppingAction =
new Geant4eSteppingAction;
289 theG4eManager->SetUserAction(theSteppingAction);
295 G4float radCyl = cDest.radius()*cm;
297 G4ThreeVector posCyl =
300 G4RotationMatrix rotCyl =
306 LogDebug(
"Geant4e") <<
"G4e - G4Rotation" << rotCyl <<
"mm";
310 G4ErrorSurfaceTarget* g4eTarget =
new G4ErrorCylSurfaceTarget(radCyl, posCyl,
314 LogDebug(
"Geant4e") <<
"G4e - Destination CMS cylinder position:" << cDest.position() <<
"cm\n"
315 <<
"G4e - Destination CMS cylinder radius:" << cDest.radius() <<
"cm\n"
316 <<
"G4e - Destination CMS cylinder rotation:" << cDest.rotation() <<
"\n";
317 LogDebug(
"Geant4e") <<
"G4e - Destination G4 cylinder position: " << posCyl <<
"mm\n"
318 <<
"G4e - Destination G4 cylinder radius:" << radCyl <<
"mm\n"
319 <<
"G4e - Destination G4 cylinder rotation:" << rotCyl <<
"\n";
327 CLHEP::Hep3Vector g4InitMom =
329 CLHEP::Hep3Vector g4InitPos =
333 LogDebug(
"Geant4e") <<
"G4e - Initial CMS point position:" << cmsInitPos
335 <<
"G4e - (Ro, eta, phi): ("
336 << cmsInitPos.
perp() <<
" cm, "
337 << cmsInitPos.
eta() <<
", "
338 << cmsInitPos.
phi().degrees() <<
" deg)\n"
339 <<
"G4e - Initial G4 point position: " << g4InitPos
340 <<
" mm, Ro = " << g4InitPos.perp() <<
" mm";
341 LogDebug(
"Geant4e") <<
"G4e - Initial CMS momentum :" << cmsInitMom
343 <<
"G4e - Initial G4 momentum : " << g4InitMom
347 int charge = ftsStart.
charge();
353 LogDebug(
"Geant4e") <<
"G4e - Particle name: " << particleName;
356 G4ErrorTrajErr g4error( 5, 1 );
361 LogDebug(
"Geant4e") <<
"G4e - Error matrix: " << g4error;
363 G4ErrorFreeTrajState* g4eTrajState =
364 new G4ErrorFreeTrajState(particleName, g4InitPos, g4InitMom, g4error);
365 LogDebug(
"Geant4e") <<
"G4e - Traj. State: " << (*g4eTrajState);
368 G4ErrorMode mode = G4ErrorMode_PropForwards;
371 mode = G4ErrorMode_PropBackwards;
372 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'backwards\'";
374 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'forwards\'";
382 mode = G4ErrorMode_PropBackwards;
383 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'backwards\'";
385 LogDebug(
"Geant4e") <<
"G4e - Propagator mode is \'forwards\'";
394 if(mode == G4ErrorMode_PropBackwards) {
398 g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
399 ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
400 g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
402 ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
404 LogDebug(
"Geant4e") <<
"G4e - Return error from propagation: " << ierr;
407 LogDebug(
"Geant4e") <<
"G4e - Error is not 0, returning invalid trajectory";
414 HepGeom::Point3D<double> posEnd = g4eTrajState->GetPosition();
415 HepGeom::Vector3D<double> momEnd = g4eTrajState->GetMomentum();
422 LogDebug(
"Geant4e") <<
"G4e - Final CMS point position:" << posEndGV
424 <<
"G4e - (Ro, eta, phi): ("
425 << posEndGV.
perp() <<
" cm, "
426 << posEndGV.
eta() <<
", "
427 << posEndGV.
phi().degrees() <<
" deg)\n"
428 <<
"G4e - Final G4 point position: " << posEnd
429 <<
" mm,\tRo =" << posEnd.perp() <<
" mm";
430 LogDebug(
"Geant4e") <<
"G4e - Final CMS momentum :" << momEndGV
432 <<
"G4e - Final G4 momentum : " << momEnd
440 G4ErrorTrajErr g4errorEnd = g4eTrajState->GetError();
443 LogDebug(
"Geant4e") <<
"G4e - Error matrix after propagation: " << g4errorEnd;
460 return propagate(ftsStart,cyl);
474 std::pair< TrajectoryStateOnSurface, double>
476 const Plane& pDest)
const {
478 theSteppingAction->reset();
483 return TsosPP(propagate(ftsStart,pDest), theSteppingAction->trackLength());
486 std::pair< TrajectoryStateOnSurface, double>
489 theSteppingAction->reset();
494 return TsosPP(propagate(ftsStart,cDest), theSteppingAction->trackLength());
497 std::pair< TrajectoryStateOnSurface, double>
499 const Plane& pDest)
const {
501 theSteppingAction->reset();
506 return TsosPP(propagate(tsosStart,pDest), theSteppingAction->trackLength());
509 std::pair< TrajectoryStateOnSurface, double>
512 theSteppingAction->reset();
517 return TsosPP(propagate(tsosStart,cDest), theSteppingAction->trackLength());
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Local3DVector LocalVector
HepGeom::Normal3D< double > globalVectorToHepNormal3D(const GlobalVector &p)
Geom::Phi< T > phi() const
float localZ(const GlobalPoint &gp) const
CLHEP::Hep3Vector globalVectorToHep3Vector(const GlobalVector &p)
TrackCharge charge() const
const CurvilinearTrajectoryError & curvilinearError() const
CLHEP::HepRotation tkRotationFToHepRotation(const TkRotation< float > &tkr)
AlgebraicSymMatrix55 g4ErrorTrajErrToAlgebraicSymMatrix55(const G4ErrorTrajErr &e, const int q)
FreeTrajectoryState const * freeState(bool withErrors=true) const
GlobalPoint hepPoint3DToGlobalPoint(const HepGeom::Point3D< double > &r)
GlobalVector momentum() const
Vector3DBase unit() const
GlobalPoint position() const
HepGeom::Point3D< double > globalPointToHepPoint3D(const GlobalPoint &r)
GlobalVector hep3VectorToGlobalVector(const CLHEP::Hep3Vector &p)
G4ErrorTrajErr algebraicSymMatrix55ToG4ErrorTrajErr(const AlgebraicSymMatrix55 &e, const int q)
CLHEP::Hep3Vector globalPointToHep3Vector(const GlobalPoint &r)