Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CascadParticle.cc
Go to the documentation of this file.
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 // $Id: G4CascadParticle.cc 67738 2013-03-05 05:54:30Z mkelsey $
27 //
28 // 20100112 M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29 // 20100114 M. Kelsey -- Replace vector<G4Double> position with G4ThreeVector
30 // 20101012 M. Kelsey -- Check for negative d2 in getPathToTheNextZone()
31 // 20110806 M. Kelsey -- Add fill() function to replicate ctor/op=() action
32 // 20110922 M. Kelsey -- Follow G4InuclParticle::print(ostream&) migration,
33 // Add stream argument to print(), add operator<<().
34 // 20111017 M. Kelsey -- Add check for zero momentum in path calculation.
35 // 20130221 M. Kelsey -- Set verbosity using global envvar, for both ctors.
36 // 20130304 M. Kelsey -- Add index data member, for use with G4CascadeHistory,
37 // and explicit copy operations
38 
39 #include "G4CascadParticle.hh"
40 #include "G4CascadeParameters.hh"
41 #include "G4ios.hh"
42 #include <cmath>
43 
44 
45 // Default constructor produces non-functional object
46 
48  : verboseLevel(G4CascadeParameters::verbose()),
49  current_zone(-1), current_path(-1.), movingIn(false),
50  reflectionCounter(0), reflected(false), generation(-1), historyId(-1) {
51  if (verboseLevel > 3) {
52  G4cout << " >>> G4CascadParticle::G4CascadParticle" << G4endl;
53  }
54 }
55 
58  const G4ThreeVector& pos, G4int izone, G4double cpath,
59  G4int gen)
60  : verboseLevel(G4CascadeParameters::verbose()),
61  theParticle(particle), position(pos),
62  current_zone(izone), current_path(cpath), movingIn(true),
63  reflectionCounter(0), reflected(false), generation(gen), historyId(-1) {
64  if (verboseLevel > 3) {
65  G4cout << " >>> G4CascadParticle::G4CascadParticle "
66  << particle.getDefinition()->GetParticleName()
67  << " @ " << pos << G4endl;
68  }
69 }
70 
71 // Copy contents, including history information
72 
74  if (&cpart != this) { // Avoid unnecessary work
75  verboseLevel = cpart.verboseLevel;
76  theParticle = cpart.theParticle;
77  position = cpart.position;
78  current_zone = cpart.current_zone;
79  current_path = cpart.current_path;
80  movingIn = cpart.movingIn;
81  reflectionCounter = cpart.reflectionCounter;
82  reflected = cpart.reflected;
83  generation = cpart.generation;
84  historyId = cpart.historyId;
85  }
86 
87  return *this;
88 }
89 
90 // Analogue to operator=() to support filling vectors w/o temporaries
91 
93  const G4ThreeVector& pos, G4int izone,
94  G4double cpath, G4int gen) {
95  if (verboseLevel > 3) G4cout << " >>> G4CascadParticle::fill" << G4endl;
96 
97  theParticle = particle;
98  position = pos;
99  current_zone = izone;
100  current_path = cpath;
101  movingIn = true;
102  reflectionCounter = 0;
103  reflected = false;
104  generation = gen;
105  historyId = -1; // New particle, no history entry yet
106 }
107 
108 
110  G4double rz_out) {
111  if (verboseLevel > 3) {
112  G4cout << " >>> G4CascadParticle::getPathToTheNextZone rz_in " << rz_in
113  << " rz_out " << rz_out << G4endl;
114  }
115 
116  const G4LorentzVector& mom = getMomentum();
117 
118  G4double path = -1.0;
119  G4double rp = mom.vect().dot(position);
120  G4double rr = position.mag2();
121  G4double pp = mom.vect().mag2();
122 
123  if (std::abs(pp) < 1e-9) { // Cut-off for "at rest" is 1 eV momentum
124  if (verboseLevel > 3) G4cout << " at rest; path length is zero" << G4endl;
125 
126  if (current_zone == 0) movingIn = false; // Allow 'stuck' to escape
127  return 0.;
128  }
129 
130  G4double ra = rr - rp * rp / pp;
131  pp = std::sqrt(pp);
132  G4double ds;
133  G4double d2;
134 
135  if (verboseLevel > 3) {
136  G4cout << " current_zone " << current_zone << " rr " << rr
137  << " rp " << rp << " pp " << pp << " ra " << ra << G4endl;
138  }
139 
140  if (current_zone == 0 || rp > 0.0) {
141  d2 = rz_out * rz_out - ra;
142  if (d2 > 0.0) {
143  ds = 1.0;
144  movingIn = false;
145  } else {
146  d2 = rz_in * rz_in - ra;
147  ds = -1.0;
148  movingIn = true;
149  }
150  } else {
151  d2 = rz_in * rz_in - ra;
152  if (d2 > 0.0) {
153  ds = -1.0;
154  movingIn = true;
155  } else {
156  d2 = rz_out * rz_out - ra;
157  ds = 1.0;
158  movingIn = false;
159  }
160  }
161 
162  if (verboseLevel > 3) G4cout << " ds " << ds << " d2 " << d2 << G4endl;
163 
164  if (d2 < 0.0 && d2 > -1e-6) d2 = 0.; // Account for round-off
165 
166  if (d2 > 0.0) path = ds * std::sqrt(d2) - rp / pp; // Avoid FPE failure
167 
168  return path;
169 }
170 
172  if (verboseLevel > 3) {
173  G4cout << " >>> G4CascadParticle::propagateAlongThePath" << G4endl;
174  }
175 
176  position += getMomentum().vect().unit()*path;
177 }
178 
179 
180 // Proper stream output (just calls print())
181 
182 std::ostream& operator<<(std::ostream& os, const G4CascadParticle& part) {
183  part.print(os);
184  return os;
185 }
186 
187 void G4CascadParticle::print(std::ostream& os) const {
188  os << " pos " << position << " zone " << current_zone
189  << " current_path " << current_path
190  << " reflectionCounter " << reflectionCounter << G4endl
191  << theParticle << G4endl;
192 }
193 
void print(std::ostream &os) const
G4LorentzVector getMomentum() const
double dot(const Hep3Vector &) const
static const G4double d2
void fill(const G4InuclElementaryParticle &particle, const G4ThreeVector &pos, G4int izone, G4double cpath, G4int gen)
const G4ParticleDefinition * getDefinition() const
void propagateAlongThePath(G4double path)
G4double getPathToTheNextZone(G4double rz_in, G4double rz_out)
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
Hep3Vector unit() const
std::ostream & operator<<(std::ostream &, const BasicVector3D< float > &)
double mag2() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4CascadParticle & operator=(const G4CascadParticle &cpart)
static const G4double pos