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 
44 Geant4ePropagator::~Geant4ePropagator() {
45 }
46 
47 //
49 //
50 
56 Geant4ePropagator::propagate (const FreeTrajectoryState& ftsStart,
57  const Plane& pDest) const {
58 
59  if(theG4eManager->PrintG4ErrorState() == "G4ErrorState_PreInit")
60  theG4eManager->InitGeant4e();
61 
62  if (!theSteppingAction) {
63  theSteppingAction = new Geant4eSteppingAction;
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 
176  if (propagationDirection() == oppositeToMomentum) {
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
272 Geant4ePropagator::propagate (const TrajectoryStateOnSurface& tsos, const Plane& plane) const {
273  const FreeTrajectoryState ftsStart = *tsos.freeState();
274  return propagate(ftsStart,plane);
275 }
276 
277 
282 Geant4ePropagator::propagate (const FreeTrajectoryState& ftsStart,
283  const Cylinder& cDest) const {
284 
285  if(theG4eManager->PrintG4ErrorState() == "G4ErrorState_PreInit")
286  theG4eManager->InitGeant4e();
287  if (!theSteppingAction) {
288  theSteppingAction = new Geant4eSteppingAction;
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 
370  if (propagationDirection() == oppositeToMomentum) {
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  LocalPoint lpos = cDest.toLocal(cmsInitPos);
380  Surface::Side theSide = cDest.side(lpos,0);
381  if(theSide==SurfaceOrientation::positiveSide){ //outside cylinder
382  mode = G4ErrorMode_PropBackwards;
383  LogDebug("Geant4e") << "G4e - Propagator mode is \'backwards\'";
384  } else { //inside cylinder
385  LogDebug("Geant4e") << "G4e - Propagator mode is \'forwards\'";
386  }
387 
388  }
389 
391  // Propagate
392 
393  int ierr;
394  if(mode == G4ErrorMode_PropBackwards) {
395  //To make geant transport the particle correctly need to give it the opposite momentum
396  //because geant flips the B field bending and adds energy instead of subtracting it
397  //but still wants the momentum "backwards"
398  g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
399  ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
400  g4eTrajState->SetMomentum( -g4eTrajState->GetMomentum());
401  } else {
402  ierr = theG4eManager->Propagate( g4eTrajState, g4eTarget, mode);
403  }
404  LogDebug("Geant4e") << "G4e - Return error from propagation: " << ierr;
405 
406  if(ierr!=0) {
407  LogDebug("Geant4e") << "G4e - Error is not 0, returning invalid trajectory";
408  return TrajectoryStateOnSurface();
409  }
410 
411  // Retrieve the state in the end from Geant4e, converte them to CMS vectors
412  // and points, and build global trajectory parameters
413  // CMS uses cm and GeV while Geant4 uses mm and MeV
414  HepGeom::Point3D<double> posEnd = g4eTrajState->GetPosition();
415  HepGeom::Vector3D<double> momEnd = g4eTrajState->GetMomentum();
416 
419 
420 
421  //DEBUG
422  LogDebug("Geant4e") << "G4e - Final CMS point position:" << posEndGV
423  << "cm\n"
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
431  << "GeV\n"
432  << "G4e - Final G4 momentum : " << momEnd
433  << " MeV";
434 
435  GlobalTrajectoryParameters tParsDest(posEndGV, momEndGV, charge, theField);
436 
437 
438  // Get the error covariance matrix from Geant4e. It comes in curvilinear
439  // coordinates so use the appropiate CMS class
440  G4ErrorTrajErr g4errorEnd = g4eTrajState->GetError();
442  curvError(TrackPropagation::g4ErrorTrajErrToAlgebraicSymMatrix55(g4errorEnd, charge));
443  LogDebug("Geant4e") << "G4e - Error matrix after propagation: " << g4errorEnd;
444 
446  // We set the SurfaceSide to atCenterOfSurface. //
448 
450 
451  return TrajectoryStateOnSurface(tParsDest, curvError, cDest, side);
452 }
453 
454 
455 //Require method with input TrajectoryStateOnSurface to be used in track fitting
456 //Don't need extra info about starting surface; use regular propagation method
458 Geant4ePropagator::propagate (const TrajectoryStateOnSurface& tsos, const Cylinder& cyl) const {
459  const FreeTrajectoryState ftsStart = *tsos.freeState();
460  return propagate(ftsStart,cyl);
461 }
462 
463 
464 //
466 //
467 
474 std::pair< TrajectoryStateOnSurface, double>
475 Geant4ePropagator::propagateWithPath (const FreeTrajectoryState& ftsStart,
476  const Plane& pDest) const {
477 
478  theSteppingAction->reset();
479 
480  //Finally build the pair<...> that needs to be returned where the second
481  //parameter is the exact path length. Currently calculated with a stepping
482  //action that adds up the length of every step
483  return TsosPP(propagate(ftsStart,pDest), theSteppingAction->trackLength());
484 }
485 
486 std::pair< TrajectoryStateOnSurface, double>
487 Geant4ePropagator::propagateWithPath (const FreeTrajectoryState& ftsStart,
488  const Cylinder& cDest) const {
489  theSteppingAction->reset();
490 
491  //Finally build the pair<...> that needs to be returned where the second
492  //parameter is the exact path length. Currently calculated with a stepping
493  //action that adds up the length of every step
494  return TsosPP(propagate(ftsStart,cDest), theSteppingAction->trackLength());
495 }
496 
497 std::pair< TrajectoryStateOnSurface, double>
498 Geant4ePropagator::propagateWithPath (const TrajectoryStateOnSurface& tsosStart,
499  const Plane& pDest) const {
500 
501  theSteppingAction->reset();
502 
503  //Finally build the pair<...> that needs to be returned where the second
504  //parameter is the exact path length. Currently calculated with a stepping
505  //action that adds up the length of every step
506  return TsosPP(propagate(tsosStart,pDest), theSteppingAction->trackLength());
507 }
508 
509 std::pair< TrajectoryStateOnSurface, double>
510 Geant4ePropagator::propagateWithPath (const TrajectoryStateOnSurface& tsosStart,
511  const Cylinder& cDest) const {
512  theSteppingAction->reset();
513 
514  //Finally build the pair<...> that needs to be returned where the second
515  //parameter is the exact path length. Currently calculated with a stepping
516  //action that adds up the length of every step
517  return TsosPP(propagate(tsosStart,cDest), theSteppingAction->trackLength());
518 }
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
Local3DVector LocalVector
Definition: LocalVector.h:12
T perp() const
Definition: PV3DBase.h:72
HepGeom::Normal3D< double > globalVectorToHepNormal3D(const GlobalVector &p)
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
float localZ(const GlobalPoint &gp) const
Definition: Plane.h:49
CLHEP::Hep3Vector globalVectorToHep3Vector(const GlobalVector &p)
PropagationDirection
TrackCharge charge() const
double charge(const std::vector< uint8_t > &Ampls)
const CurvilinearTrajectoryError & curvilinearError() const
Definition: Plane.h:17
tuple field
Definition: statics.py:62
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
Definition: Vector3DBase.h:57
GlobalPoint position() const
HepGeom::Point3D< double > globalPointToHepPoint3D(const GlobalPoint &r)
GlobalVector hep3VectorToGlobalVector(const CLHEP::Hep3Vector &p)
T eta() const
Definition: PV3DBase.h:76
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
G4ErrorTrajErr algebraicSymMatrix55ToG4ErrorTrajErr(const AlgebraicSymMatrix55 &e, const int q)
tuple cout
Definition: gather_cfg.py:121
dbl *** dir
Definition: mlp_gen.cc:35
CLHEP::Hep3Vector globalPointToHep3Vector(const GlobalPoint &r)