CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Geant4ePropagator.cc
Go to the documentation of this file.
1 
2 //Geant4e
6 
7 //CMSSW
15 
16 //Geant4
17 #include "G4ErrorFreeTrajState.hh"
18 #include "G4ErrorPlaneSurfaceTarget.hh"
19 #include "G4ErrorCylSurfaceTarget.hh"
20 #include "G4ErrorPropagatorData.hh"
21 #include "G4EventManager.hh"
22 #include "G4SteppingControl.hh"
23 
24 //CLHEP
25 #include "CLHEP/Units/GlobalSystemOfUnits.h"
26 
27 
31  const char* particleName,
33  Propagator(dir),
34  theField(field),
35  theParticleName(particleName),
36  theG4eManager(G4ErrorPropagatorManager::GetErrorPropagatorManager()),
37  theSteppingAction(0) {
38 
39  G4ErrorPropagatorData::SetVerbose(0);
40 }
41 
45 }
46 
47 //
49 //
50 
57  const Plane& pDest) const {
58 
59  if(theG4eManager->PrintG4ErrorState() == "G4ErrorState_PreInit")
60  theG4eManager->InitGeant4e();
61 
62  if (!theSteppingAction) {
64  theG4eManager->SetUserAction(theSteppingAction);
65  }
66 
68  // Construct the target surface
69  //
70 
71  //* Get position and normal (orientation) of the destination plane
72  GlobalPoint posPlane = pDest.toGlobal(LocalPoint(0,0,0));
73  GlobalVector normalPlane = pDest.toGlobal(LocalVector(0,0,1.));
74  normalPlane = normalPlane.unit();
75 
76  //* Transform this into HepGeom::Point3D<double> and HepGeom::Normal3D<double> that define a plane for
77  // Geant4e.
78  // CMS uses cm and GeV while Geant4 uses mm and MeV
79  HepGeom::Point3D<double> surfPos =
81  HepGeom::Normal3D<double> surfNorm =
83 
84  //DEBUG
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 : "
95  << normalPlane;
96  LogDebug("Geant4e") << "G4e - Distance from plane position to plane: "
97  << pDest.localZ(posPlane) << " cm";
98  //DEBUG
99 
100  //* Set the target surface
101  G4ErrorSurfaceTarget* g4eTarget = new G4ErrorPlaneSurfaceTarget(surfNorm,
102  surfPos);
103 
104  //g4eTarget->Dump("G4e - ");
105  //
107 
109  // Find initial point
110  //
111 
112  // * Get the starting point and direction and convert them to CLHEP::Hep3Vector
113  // for G4. CMS uses cm and GeV while Geant4 uses mm and MeV
114  GlobalPoint cmsInitPos = ftsStart.position();
115  GlobalVector cmsInitMom = ftsStart.momentum();
116 
117  CLHEP::Hep3Vector g4InitPos =
119  CLHEP::Hep3Vector g4InitMom =
121 
122  //DEBUG
123  LogDebug("Geant4e") << "G4e - Initial CMS point position:" << cmsInitPos
124  << "cm\n"
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
132  << "GeV\n"
133  << "G4e - Initial G4 momentum : " << g4InitMom
134  << " MeV";
135  LogDebug("Geant4e") << "G4e - Distance from initial point to plane: "
136  << pDest.localZ(cmsInitPos) << " cm";
137  //DEBUG
138 
139  //
141 
143  // Set particle name
144  //
145  int charge = ftsStart.charge();
146  std::string particleName = theParticleName;
147 
148  if (charge > 0) {
149  particleName += "+";
150  } else {
151  particleName += "-";
152  }
153 
154  LogDebug("Geant4e") << "G4e - Particle name: " << particleName;
155 
156  //
158 
160  //Set the error and trajectories, and finally propagate
161  //
162  G4ErrorTrajErr g4error( 5, 1 );
163  if(ftsStart.hasError()) {
164  const CurvilinearTrajectoryError initErr = ftsStart.curvilinearError();
165  g4error = TrackPropagation::algebraicSymMatrix55ToG4ErrorTrajErr( initErr , charge); //The error matrix
166  }
167  LogDebug("Geant4e") << "G4e - Error matrix: " << g4error;
168 
169  G4ErrorFreeTrajState* g4eTrajState =
170  new G4ErrorFreeTrajState(particleName, g4InitPos, g4InitMom, g4error);
171  LogDebug("Geant4e") << "G4e - Traj. State: " << (*g4eTrajState);
172 
173  //Set the mode of propagation according to the propagation direction
174  G4ErrorMode mode = G4ErrorMode_PropForwards;
175 
177  mode = G4ErrorMode_PropBackwards;
178  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\'";
179  } else if(propagationDirection() == alongMomentum) {
180  LogDebug("Geant4e") << "G4e - Propagator mode is \'forwards\'";
181  } else { //Mode must be anyDirection then - need to figure out for Geant which it is
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;
186  } else {
187  mode = G4ErrorMode_PropBackwards;
188  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\'";
189  std::cout << ", got backwards" << std::endl;
190  }
191  }
192 
193  //
195 
197  // Propagate
198 
199  int ierr;
200  if(mode == G4ErrorMode_PropBackwards) {
201  //To make geant transport the particle correctly need to give it the opposite momentum
202  //because geant flips the B field bending and adds energy instead of subtracting it
203  //but still wants the momentum "backwards"
204  g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
205  ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
206  g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
207  } else {
208  ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
209  }
210  LogDebug("Geant4e") << "G4e - Return error from propagation: " << ierr;
211 
212  if(ierr!=0) {
213  LogDebug("Geant4e") << "G4e - Error is not 0, returning invalid trajectory";
214  return TrajectoryStateOnSurface();
215  }
216 
217  //
219 
221  // Retrieve the state in the end from Geant4e, convert them to CMS vectors
222  // and points, and build global trajectory parameters.
223  // CMS uses cm and GeV while Geant4 uses mm and MeV
224  //
225  HepGeom::Point3D<double> posEnd = g4eTrajState->GetPosition();
226  HepGeom::Vector3D<double> momEnd = g4eTrajState->GetMomentum();
227 
230 
231  //DEBUG
232  LogDebug("Geant4e") << "G4e - Final CMS point position:" << posEndGV
233  << "cm\n"
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
241  << "GeV\n"
242  << "G4e - Final G4 momentum : " << momEnd
243  << " MeV";
244  LogDebug("Geant4e") << "G4e - Distance from final point to plane: "
245  << pDest.localZ(posEndGV) << " cm";
246  //DEBUG
247 
248  GlobalTrajectoryParameters tParsDest(posEndGV, momEndGV, charge, theField);
249 
250 
251  // Get the error covariance matrix from Geant4e. It comes in curvilinear
252  // coordinates so use the appropiate CMS class
253  G4ErrorTrajErr g4errorEnd = g4eTrajState->GetError();
255  curvError(TrackPropagation::g4ErrorTrajErrToAlgebraicSymMatrix55(g4errorEnd, charge));
256  LogDebug("Geant4e") << "G4e - Error matrix after propagation: " << g4errorEnd;
257 
259  // We set the SurfaceSide to atCenterOfSurface. //
261  LogDebug("Geant4e") << "G4e - SurfaceSide is always atCenterOfSurface after propagation";
263  //
265 
266  return TrajectoryStateOnSurface(tParsDest, curvError, pDest, side);
267 }
268 
269 //Require method with input TrajectoryStateOnSurface to be used in track fitting
270 //Don't need extra info about starting surface; use regular propagation method
273  const FreeTrajectoryState ftsStart = *tsos.freeState();
274  return propagate(ftsStart,plane);
275 }
276 
277 
283  const Cylinder& cDest) const {
284 
285  if(theG4eManager->PrintG4ErrorState() == "G4ErrorState_PreInit")
286  theG4eManager->InitGeant4e();
287  if (!theSteppingAction) {
289  theG4eManager->SetUserAction(theSteppingAction);
290  }
291 
292  //Get Cylinder parameters.
293  //CMS uses cm and GeV while Geant4 uses mm and MeV.
294  // - Radius
295  G4float radCyl = cDest.radius()*cm;
296  // - Position: PositionType & GlobalPoint are Basic3DPoint<float,GlobalTag>
297  G4ThreeVector posCyl =
299  // - Rotation: Type in CMSSW is RotationType == TkRotation<T>, T=float
300  G4RotationMatrix rotCyl =
302 
303  //DEBUG
304  TkRotation<float> rotation = cDest.rotation();
305  LogDebug("Geant4e") << "G4e - TkRotation" << rotation;
306  LogDebug("Geant4e") << "G4e - G4Rotation" << rotCyl << "mm";
307 
308 
309  //Set the target surface
310  G4ErrorSurfaceTarget* g4eTarget = new G4ErrorCylSurfaceTarget(radCyl, posCyl,
311  rotCyl);
312 
313  //DEBUG
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";
320 
321 
322  //Get the starting point and direction and convert them to CLHEP::Hep3Vector for G4
323  //CMS uses cm and GeV while Geant4 uses mm and MeV
324  GlobalPoint cmsInitPos = ftsStart.position();
325  GlobalVector cmsInitMom = ftsStart.momentum();
326 
327  CLHEP::Hep3Vector g4InitMom =
329  CLHEP::Hep3Vector g4InitPos =
331 
332  //DEBUG
333  LogDebug("Geant4e") << "G4e - Initial CMS point position:" << cmsInitPos
334  << "cm\n"
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
342  << "GeV\n"
343  << "G4e - Initial G4 momentum : " << g4InitMom
344  << " MeV";
345 
346  //Set particle name
347  int charge = ftsStart.charge();
348  std::string particleName = theParticleName;
349  if (charge > 0)
350  particleName += "+";
351  else
352  particleName += "-";
353  LogDebug("Geant4e") << "G4e - Particle name: " << particleName;
354 
355  //Set the error and trajectories, and finally propagate
356  G4ErrorTrajErr g4error( 5, 1 );
357  if(ftsStart.hasError()) {
358  const CurvilinearTrajectoryError initErr = ftsStart.curvilinearError();
359  g4error = TrackPropagation::algebraicSymMatrix55ToG4ErrorTrajErr( initErr , charge); //The error matrix
360  }
361  LogDebug("Geant4e") << "G4e - Error matrix: " << g4error;
362 
363  G4ErrorFreeTrajState* g4eTrajState =
364  new G4ErrorFreeTrajState(particleName, g4InitPos, g4InitMom, g4error);
365  LogDebug("Geant4e") << "G4e - Traj. State: " << (*g4eTrajState);
366 
367  //Set the mode of propagation according to the propagation direction
368  G4ErrorMode mode = G4ErrorMode_PropForwards;
369 
371  mode = G4ErrorMode_PropBackwards;
372  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\'";
373  } else if(propagationDirection() == alongMomentum) {
374  LogDebug("Geant4e") << "G4e - Propagator mode is \'forwards\'";
375  } else {
376  //------------------------------------
377  //For cylinder assume outside is backwards, inside is along
378  //General use for particles from collisions
379  LocalVector lmom = cDest.toLocal(cmsInitMom);
380  LocalPoint lpos = cDest.toLocal(cmsInitPos);
381  Surface::Side theSide = cDest.side(lpos,0);
382  if(theSide==SurfaceOrientation::positiveSide){ //outside cylinder
383  mode = G4ErrorMode_PropBackwards;
384  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\'";
385  } else { //inside cylinder
386  LogDebug("Geant4e") << "G4e - Propagator mode is \'forwards\'";
387  }
388 
389  }
390 
392  // Propagate
393 
394  int ierr;
395  if(mode == G4ErrorMode_PropBackwards) {
396  //To make geant transport the particle correctly need to give it the opposite momentum
397  //because geant flips the B field bending and adds energy instead of subtracting it
398  //but still wants the momentum "backwards"
399  g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
400  ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
401  g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
402  } else {
403  ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
404  }
405  LogDebug("Geant4e") << "G4e - Return error from propagation: " << ierr;
406 
407  if(ierr!=0) {
408  LogDebug("Geant4e") << "G4e - Error is not 0, returning invalid trajectory";
409  return TrajectoryStateOnSurface();
410  }
411 
412  // Retrieve the state in the end from Geant4e, converte them to CMS vectors
413  // and points, and build global trajectory parameters
414  // CMS uses cm and GeV while Geant4 uses mm and MeV
415  HepGeom::Point3D<double> posEnd = g4eTrajState->GetPosition();
416  HepGeom::Vector3D<double> momEnd = g4eTrajState->GetMomentum();
417 
420 
421 
422  //DEBUG
423  LogDebug("Geant4e") << "G4e - Final CMS point position:" << posEndGV
424  << "cm\n"
425  << "G4e - (Ro, eta, phi): ("
426  << posEndGV.perp() << " cm, "
427  << posEndGV.eta() << ", "
428  << posEndGV.phi().degrees() << " deg)\n"
429  << "G4e - Final G4 point position: " << posEnd
430  << " mm,\tRo =" << posEnd.perp() << " mm";
431  LogDebug("Geant4e") << "G4e - Final CMS momentum :" << momEndGV
432  << "GeV\n"
433  << "G4e - Final G4 momentum : " << momEnd
434  << " MeV";
435 
436  GlobalTrajectoryParameters tParsDest(posEndGV, momEndGV, charge, theField);
437 
438 
439  // Get the error covariance matrix from Geant4e. It comes in curvilinear
440  // coordinates so use the appropiate CMS class
441  G4ErrorTrajErr g4errorEnd = g4eTrajState->GetError();
443  curvError(TrackPropagation::g4ErrorTrajErrToAlgebraicSymMatrix55(g4errorEnd, charge));
444  LogDebug("Geant4e") << "G4e - Error matrix after propagation: " << g4errorEnd;
445 
447  // We set the SurfaceSide to atCenterOfSurface. //
449 
451 
452  return TrajectoryStateOnSurface(tParsDest, curvError, cDest, side);
453 }
454 
455 
456 //Require method with input TrajectoryStateOnSurface to be used in track fitting
457 //Don't need extra info about starting surface; use regular propagation method
460  const FreeTrajectoryState ftsStart = *tsos.freeState();
461  return propagate(ftsStart,cyl);
462 }
463 
464 
465 //
467 //
468 
475 std::pair< TrajectoryStateOnSurface, double>
477  const Plane& pDest) const {
478 
480 
481  //Finally build the pair<...> that needs to be returned where the second
482  //parameter is the exact path length. Currently calculated with a stepping
483  //action that adds up the length of every step
484  return TsosPP(propagate(ftsStart,pDest), theSteppingAction->trackLength());
485 }
486 
487 std::pair< TrajectoryStateOnSurface, double>
489  const Cylinder& cDest) const {
491 
492  //Finally build the pair<...> that needs to be returned where the second
493  //parameter is the exact path length. Currently calculated with a stepping
494  //action that adds up the length of every step
495  return TsosPP(propagate(ftsStart,cDest), theSteppingAction->trackLength());
496 }
497 
498 std::pair< TrajectoryStateOnSurface, double>
500  const Plane& pDest) const {
501 
503 
504  //Finally build the pair<...> that needs to be returned where the second
505  //parameter is the exact path length. Currently calculated with a stepping
506  //action that adds up the length of every step
507  return TsosPP(propagate(tsosStart,pDest), theSteppingAction->trackLength());
508 }
509 
510 std::pair< TrajectoryStateOnSurface, double>
512  const Cylinder& cDest) const {
514 
515  //Finally build the pair<...> that needs to be returned where the second
516  //parameter is the exact path length. Currently calculated with a stepping
517  //action that adds up the length of every step
518  return TsosPP(propagate(tsosStart,cDest), theSteppingAction->trackLength());
519 }
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:78
std::pair< TrajectoryStateOnSurface, double > TsosPP
Local3DVector LocalVector
Definition: LocalVector.h:12
T perp() const
Definition: PV3DBase.h:66
virtual PropagationDirection propagationDirection() const
Definition: Propagator.h:143
HepGeom::Normal3D< double > globalVectorToHepNormal3D(const GlobalVector &p)
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
virtual std::pair< TrajectoryStateOnSurface, double > propagateWithPath(const FreeTrajectoryState &, const Plane &) const
float localZ(const GlobalPoint &gp) const
Fast access to distance from plane for a point.
Definition: Plane.h:52
std::string theParticleName
CLHEP::Hep3Vector globalVectorToHep3Vector(const GlobalVector &p)
PropagationDirection
TrackCharge charge() const
const MagneticField * theField
double charge(const std::vector< uint8_t > &Ampls)
const CurvilinearTrajectoryError & curvilinearError() const
Definition: Plane.h:17
G4ErrorPropagatorManager * theG4eManager
CLHEP::HepRotation tkRotationFToHepRotation(const TkRotation< float > &tkr)
virtual Side side(const LocalPoint &p, Scalar toler) const
Definition: Cylinder.cc:9
Geant4ePropagator(const MagneticField *field=0, const char *particleName="mu", PropagationDirection dir=alongMomentum)
Scalar radius() const
Radius of the cylinder.
Definition: Cylinder.h:55
AlgebraicSymMatrix55 g4ErrorTrajErrToAlgebraicSymMatrix55(const G4ErrorTrajErr &e, const int q)
FreeTrajectoryState * freeState(bool withErrors=true) const
LocalPoint toLocal(const GlobalPoint &gp) const
GlobalPoint hepPoint3DToGlobalPoint(const HepGeom::Point3D< double > &r)
GlobalVector momentum() const
Geant4eSteppingAction * theSteppingAction
Vector3DBase unit() const
Definition: Vector3DBase.h:57
GlobalPoint position() const
HepGeom::Point3D< double > globalPointToHepPoint3D(const GlobalPoint &r)
int mode
Definition: AMPTWrapper.h:139
virtual ~Geant4ePropagator()
GlobalVector hep3VectorToGlobalVector(const CLHEP::Hep3Vector &p)
T eta() const
Definition: PV3DBase.h:70
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
const RotationType & rotation() const
G4ErrorTrajErr algebraicSymMatrix55ToG4ErrorTrajErr(const AlgebraicSymMatrix55 &e, const int q)
tuple cout
Definition: gather_cfg.py:41
virtual TrajectoryStateOnSurface propagate(const FreeTrajectoryState &ftsStart, const Plane &pDest) const
dbl *** dir
Definition: mlp_gen.cc:35
T degrees() const
Definition: Phi.h:51
const PositionType & position() const
CLHEP::Hep3Vector globalPointToHep3Vector(const GlobalPoint &r)