| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /** \file G4INCLIntersection.hh
39 * \brief Simple class for computing intersections between a straight line and a sphere.
40 *
41 * \date 12 December 2011
42 * \author Davide Mancusi
43 */
44
45 #ifndef G4INCLINTERSECTION_HH
46 #define G4INCLINTERSECTION_HH 1
47
48 #include "G4INCLThreeVector.hh"
49 #include <utility>
50
51 namespace G4INCL {
52
53 /** \brief Intersection-point structure
54 *
55 * The structure contains the time and position of the intersection point
56 * of a trajectory with a surface, if it exists.
57 */
58 struct Intersection {
59 Intersection(const G4bool e, const G4double t, const ThreeVector &p) :
60 exists(e), time(t), position(p) {}
61 G4bool exists;
62 G4double time;
63 ThreeVector position;
64 };
65
66 namespace IntersectionFactory {
67
68 /** \brief Compute the first intersection of a straight particle
69 * trajectory with a sphere.
70 *
71 * \param x0 the starting position of the trajectory
72 * \param p the trajectory direction
73 * \param r the radius of the sphere (centred in the origin)
74 * \return an Intersection. The G4bool is true if an intersection exists,
75 * in which case its position is stored in the ThreeVector and
76 * its time in the G4double.
77 */
78 Intersection getEarlierTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r);
79 /** \brief Compute the second intersection of a straight particle
80 * trajectory with a sphere.
81 *
82 * \param x0 the starting position of the trajectory
83 * \param p the trajectory direction
84 * \param r the radius of the sphere (centred in the origin)
85 * \return an Intersection. The G4bool is true if an intersection exists,
86 * in which case its position is stored in the ThreeVector and
87 * its time in the G4double.
88 */
89 Intersection getLaterTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r);
90 /** \brief Compute both intersections of a straight particle
91 * trajectory with a sphere.
92 *
93 * \param x0 the starting position of the trajectory
94 * \param p the trajectory direction
95 * \param r the radius of the sphere (centred in the origin)
96 * \return an Intersection. The G4bool is true if an intersection exists,
97 * in which case its position is stored in the ThreeVector and
98 * its time in the G4double.
99 */
100 std::pair<Intersection,Intersection> getTrajectoryIntersections(const ThreeVector &x0, const ThreeVector &p, const G4double r);
101
102 namespace {
103
104 Intersection getTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &v, const G4double r, const G4bool earliest) {
105 const G4double scalarVelocity = v.mag();
106 ThreeVector velocityUnitVector = v / scalarVelocity;
107
108 ThreeVector positionTransverse = x0 - velocityUnitVector * x0.dot(velocityUnitVector);
109 const G4double impactParameter = positionTransverse.mag();
110
111 const G4double r2 = r*r;
112 G4double distanceZ2 = r2 - impactParameter * impactParameter;
113 if(distanceZ2 < 0.0)
114 return Intersection(false, 0.0, ThreeVector());
115
116 const G4double distanceZ = std::sqrt(distanceZ2);
117 const ThreeVector position = positionTransverse + velocityUnitVector * (earliest ? -distanceZ : distanceZ);
118 const G4double time = (position-x0).dot(velocityUnitVector)/scalarVelocity;
119 return Intersection(true, time, position);
120 }
121
122 }
123
124 inline Intersection getEarlierTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r) {
125 return getTrajectoryIntersection(x0, p, r, true);
126 }
127
128 inline Intersection getLaterTrajectoryIntersection(const ThreeVector &x0, const ThreeVector &p, const G4double r) {
129 return getTrajectoryIntersection(x0, p, r, false);
130 }
131
132 inline std::pair<Intersection,Intersection> getTrajectoryIntersections(const ThreeVector &x0, const ThreeVector &p, const G4double r) {
133 return std::make_pair(
134 getTrajectoryIntersection(x0, p, r, true),
135 getTrajectoryIntersection(x0, p, r, false)
136 );
137 }
138 }
139
140 }
141
142 #endif /* G4INCLINTERSECTION_HH */
143