| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // GEANT 4 class implementation file 28 // GEANT 4 class implementation file
29 // 29 //
30 // CERN, Geneva, Switzerland 30 // CERN, Geneva, Switzerland
31 // 31 //
32 // File name: G4RKPropagation.cc 32 // File name: G4RKPropagation.cc
33 // 33 //
34 // Author: Alessandro Brunengo (Al 34 // Author: Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
35 // 35 //
36 // Creation date: 6 June 2000 36 // Creation date: 6 June 2000
37 // ------------------------------------------- 37 // -------------------------------------------------------------------
38 #include "G4RKPropagation.hh" 38 #include "G4RKPropagation.hh"
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 // nuclear fields 41 // nuclear fields
42 #include "G4VNuclearField.hh" 42 #include "G4VNuclearField.hh"
43 #include "G4ProtonField.hh" 43 #include "G4ProtonField.hh"
44 #include "G4NeutronField.hh" 44 #include "G4NeutronField.hh"
45 #include "G4AntiProtonField.hh" 45 #include "G4AntiProtonField.hh"
46 #include "G4KaonPlusField.hh" 46 #include "G4KaonPlusField.hh"
47 #include "G4KaonMinusField.hh" 47 #include "G4KaonMinusField.hh"
48 #include "G4KaonZeroField.hh" 48 #include "G4KaonZeroField.hh"
49 #include "G4PionPlusField.hh" 49 #include "G4PionPlusField.hh"
50 #include "G4PionMinusField.hh" 50 #include "G4PionMinusField.hh"
51 #include "G4PionZeroField.hh" 51 #include "G4PionZeroField.hh"
52 #include "G4SigmaPlusField.hh" 52 #include "G4SigmaPlusField.hh"
53 #include "G4SigmaMinusField.hh" 53 #include "G4SigmaMinusField.hh"
54 #include "G4SigmaZeroField.hh" 54 #include "G4SigmaZeroField.hh"
55 // particles properties 55 // particles properties
56 #include "G4Proton.hh" 56 #include "G4Proton.hh"
57 #include "G4Neutron.hh" 57 #include "G4Neutron.hh"
58 #include "G4AntiProton.hh" 58 #include "G4AntiProton.hh"
59 #include "G4KaonPlus.hh" 59 #include "G4KaonPlus.hh"
60 #include "G4KaonMinus.hh" 60 #include "G4KaonMinus.hh"
61 #include "G4KaonZero.hh" 61 #include "G4KaonZero.hh"
62 #include "G4PionPlus.hh" 62 #include "G4PionPlus.hh"
63 #include "G4PionMinus.hh" 63 #include "G4PionMinus.hh"
64 #include "G4PionZero.hh" 64 #include "G4PionZero.hh"
65 #include "G4SigmaPlus.hh" 65 #include "G4SigmaPlus.hh"
66 #include "G4SigmaMinus.hh" 66 #include "G4SigmaMinus.hh"
67 #include "G4SigmaZero.hh" 67 #include "G4SigmaZero.hh"
68 68
69 #include "globals.hh" 69 #include "globals.hh"
70 70
71 #include "G4KM_OpticalEqRhs.hh" 71 #include "G4KM_OpticalEqRhs.hh"
72 #include "G4KM_NucleonEqRhs.hh" 72 #include "G4KM_NucleonEqRhs.hh"
73 #include "G4ClassicalRK4.hh" 73 #include "G4ClassicalRK4.hh"
74 #include "G4MagIntegratorDriver.hh" 74 #include "G4MagIntegratorDriver.hh"
75 75
76 #include "G4LorentzRotation.hh" 76 #include "G4LorentzRotation.hh"
77 77
78 // unsigned EncodingHashFun(const G4int& aEnco 78 // unsigned EncodingHashFun(const G4int& aEncoding);
79 79
80 G4RKPropagation::G4RKPropagation() : 80 G4RKPropagation::G4RKPropagation() :
81 theOuterRadius(0), theNucleus(0), 81 theOuterRadius(0), theNucleus(0),
82 theFieldMap(0), theEquationMap(0), 82 theFieldMap(0), theEquationMap(0),
83 theField(0) 83 theField(0)
84 { } 84 { }
85 85
86 86
87 G4RKPropagation::~G4RKPropagation() 87 G4RKPropagation::~G4RKPropagation()
88 { 88 {
89 // free theFieldMap memory 89 // free theFieldMap memory
90 if(theFieldMap) delete_FieldsAndMap(theFiel 90 if(theFieldMap) delete_FieldsAndMap(theFieldMap);
91 91
92 // free theEquationMap memory 92 // free theEquationMap memory
93 if(theEquationMap) delete_EquationsAndMap(t 93 if(theEquationMap) delete_EquationsAndMap(theEquationMap);
94 94
95 if (theField) delete theField; 95 if (theField) delete theField;
96 } 96 }
97 97
98 //-------------------------------------------- 98 //----------------------------------------------------------------------------
99 void G4RKPropagation::Init(G4V3DNucleus * nucl 99 void G4RKPropagation::Init(G4V3DNucleus * nucleus)
100 //-------------------------------------------- 100 //----------------------------------------------------------------------------
101 { 101 {
102 102
103 // free theFieldMap memory 103 // free theFieldMap memory
104 if(theFieldMap) delete_FieldsAndMap(theFiel 104 if(theFieldMap) delete_FieldsAndMap(theFieldMap);
105 105
106 // free theEquationMap memory 106 // free theEquationMap memory
107 if(theEquationMap) delete_EquationsAndMap(t 107 if(theEquationMap) delete_EquationsAndMap(theEquationMap);
108 108
109 if (theField) delete theField; 109 if (theField) delete theField;
110 110
111 // Initialize the nuclear field map. 111 // Initialize the nuclear field map.
112 theNucleus = nucleus; 112 theNucleus = nucleus;
113 theOuterRadius = theNucleus->GetOuterRadius 113 theOuterRadius = theNucleus->GetOuterRadius();
114 114
115 theFieldMap = new std::map <G4int, G4VNucle 115 theFieldMap = new std::map <G4int, G4VNuclearField*, std::less<G4int> >;
116 116
117 (*theFieldMap)[G4Proton::Proton()->GetPDGEn 117 (*theFieldMap)[G4Proton::Proton()->GetPDGEncoding()] = new G4ProtonField(theNucleus);
118 (*theFieldMap)[G4Neutron::Neutron()->GetPDG 118 (*theFieldMap)[G4Neutron::Neutron()->GetPDGEncoding()] = new G4NeutronField(theNucleus);
119 (*theFieldMap)[G4AntiProton::AntiProton()-> 119 (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = new G4AntiProtonField(theNucleus);
120 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetP 120 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = new G4KaonPlusField(theNucleus);
121 (*theFieldMap)[G4KaonMinus::KaonMinus()->Ge 121 (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = new G4KaonMinusField(theNucleus);
122 (*theFieldMap)[G4KaonZero::KaonZero()->GetP 122 (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = new G4KaonZeroField(theNucleus);
123 (*theFieldMap)[G4PionPlus::PionPlus()->GetP 123 (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = new G4PionPlusField(theNucleus);
124 (*theFieldMap)[G4PionMinus::PionMinus()->Ge 124 (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = new G4PionMinusField(theNucleus);
125 (*theFieldMap)[G4PionZero::PionZero()->GetP 125 (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()] = new G4PionZeroField(theNucleus);
126 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->Ge 126 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = new G4SigmaPlusField(theNucleus);
127 (*theFieldMap)[G4SigmaMinus::SigmaMinus()-> 127 (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = new G4SigmaMinusField(theNucleus);
128 (*theFieldMap)[G4SigmaZero::SigmaZero()->Ge 128 (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = new G4SigmaZeroField(theNucleus);
129 129
130 theEquationMap = new std::map <G4int, G4Mag 130 theEquationMap = new std::map <G4int, G4Mag_EqRhs*, std::less<G4int> >;
131 131
132 // theField needed by the design of G4Mag_e 132 // theField needed by the design of G4Mag_eqRhs
133 theField = new G4KM_DummyField; //Field 133 theField = new G4KM_DummyField; //Field not needed for integration
134 G4KM_OpticalEqRhs * opticalEq; 134 G4KM_OpticalEqRhs * opticalEq;
135 G4KM_NucleonEqRhs * nucleonEq; 135 G4KM_NucleonEqRhs * nucleonEq;
136 G4double mass; 136 G4double mass;
137 G4double opticalCoeff; 137 G4double opticalCoeff;
138 138
139 nucleonEq = new G4KM_NucleonEqRhs(theField, 139 nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus);
140 mass = G4Proton::Proton()->GetPDGMass(); 140 mass = G4Proton::Proton()->GetPDGMass();
141 nucleonEq->SetMass(mass); 141 nucleonEq->SetMass(mass);
142 (*theEquationMap)[G4Proton::Proton()->GetPD 142 (*theEquationMap)[G4Proton::Proton()->GetPDGEncoding()] = nucleonEq;
143 143
144 nucleonEq = new G4KM_NucleonEqRhs(theField, 144 nucleonEq = new G4KM_NucleonEqRhs(theField, theNucleus);
145 mass = G4Neutron::Neutron()->GetPDGMass(); 145 mass = G4Neutron::Neutron()->GetPDGMass();
146 nucleonEq->SetMass(mass); 146 nucleonEq->SetMass(mass);
147 (*theEquationMap)[G4Neutron::Neutron()->Get 147 (*theEquationMap)[G4Neutron::Neutron()->GetPDGEncoding()] = nucleonEq;
148 148
149 opticalEq = new G4KM_OpticalEqRhs(theField, 149 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
150 mass = G4AntiProton::AntiProton()->GetPDGMa 150 mass = G4AntiProton::AntiProton()->GetPDGMass();
151 opticalCoeff = 151 opticalCoeff =
152 (*theFieldMap)[G4AntiProton::AntiProt 152 (*theFieldMap)[G4AntiProton::AntiProton()->GetPDGEncoding()]->GetCoeff();
153 opticalEq->SetFactor(mass,opticalCoeff); 153 opticalEq->SetFactor(mass,opticalCoeff);
154 (*theEquationMap)[G4AntiProton::AntiProton( 154 (*theEquationMap)[G4AntiProton::AntiProton()->GetPDGEncoding()] = opticalEq;
155 155
156 opticalEq = new G4KM_OpticalEqRhs(theField, 156 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
157 mass = G4KaonPlus::KaonPlus()->GetPDGMass() 157 mass = G4KaonPlus::KaonPlus()->GetPDGMass();
158 opticalCoeff = 158 opticalCoeff =
159 (*theFieldMap)[G4KaonPlus::KaonPlus() 159 (*theFieldMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()]->GetCoeff();
160 opticalEq->SetFactor(mass,opticalCoeff); 160 opticalEq->SetFactor(mass,opticalCoeff);
161 (*theEquationMap)[G4KaonPlus::KaonPlus()->G 161 (*theEquationMap)[G4KaonPlus::KaonPlus()->GetPDGEncoding()] = opticalEq;
162 162
163 opticalEq = new G4KM_OpticalEqRhs(theField, 163 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
164 mass = G4KaonMinus::KaonMinus()->GetPDGMass 164 mass = G4KaonMinus::KaonMinus()->GetPDGMass();
165 opticalCoeff = 165 opticalCoeff =
166 (*theFieldMap)[G4KaonMinus::KaonMinus 166 (*theFieldMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()]->GetCoeff();
167 opticalEq->SetFactor(mass,opticalCoeff); 167 opticalEq->SetFactor(mass,opticalCoeff);
168 (*theEquationMap)[G4KaonMinus::KaonMinus()- 168 (*theEquationMap)[G4KaonMinus::KaonMinus()->GetPDGEncoding()] = opticalEq;
169 169
170 opticalEq = new G4KM_OpticalEqRhs(theField, 170 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
171 mass = G4KaonZero::KaonZero()->GetPDGMass() 171 mass = G4KaonZero::KaonZero()->GetPDGMass();
172 opticalCoeff = 172 opticalCoeff =
173 (*theFieldMap)[G4KaonZero::KaonZero() 173 (*theFieldMap)[G4KaonZero::KaonZero()->GetPDGEncoding()]->GetCoeff();
174 opticalEq->SetFactor(mass,opticalCoeff); 174 opticalEq->SetFactor(mass,opticalCoeff);
175 (*theEquationMap)[G4KaonZero::KaonZero()->G 175 (*theEquationMap)[G4KaonZero::KaonZero()->GetPDGEncoding()] = opticalEq;
176 176
177 opticalEq = new G4KM_OpticalEqRhs(theField, 177 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
178 mass = G4PionPlus::PionPlus()->GetPDGMass() 178 mass = G4PionPlus::PionPlus()->GetPDGMass();
179 opticalCoeff = 179 opticalCoeff =
180 (*theFieldMap)[G4PionPlus::PionPlus() 180 (*theFieldMap)[G4PionPlus::PionPlus()->GetPDGEncoding()]->GetCoeff();
181 opticalEq->SetFactor(mass,opticalCoeff); 181 opticalEq->SetFactor(mass,opticalCoeff);
182 (*theEquationMap)[G4PionPlus::PionPlus()->G 182 (*theEquationMap)[G4PionPlus::PionPlus()->GetPDGEncoding()] = opticalEq;
183 183
184 opticalEq = new G4KM_OpticalEqRhs(theField, 184 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
185 mass = G4PionMinus::PionMinus()->GetPDGMass 185 mass = G4PionMinus::PionMinus()->GetPDGMass();
186 opticalCoeff = 186 opticalCoeff =
187 (*theFieldMap)[G4PionMinus::PionMinus 187 (*theFieldMap)[G4PionMinus::PionMinus()->GetPDGEncoding()]->GetCoeff();
188 opticalEq->SetFactor(mass,opticalCoeff); 188 opticalEq->SetFactor(mass,opticalCoeff);
189 (*theEquationMap)[G4PionMinus::PionMinus()- 189 (*theEquationMap)[G4PionMinus::PionMinus()->GetPDGEncoding()] = opticalEq;
190 190
191 opticalEq = new G4KM_OpticalEqRhs(theField, 191 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
192 mass = G4PionZero::PionZero()->GetPDGMass() 192 mass = G4PionZero::PionZero()->GetPDGMass();
193 opticalCoeff = 193 opticalCoeff =
194 (*theFieldMap)[G4PionZero::PionZero() 194 (*theFieldMap)[G4PionZero::PionZero()->GetPDGEncoding()]->GetCoeff();
195 opticalEq->SetFactor(mass,opticalCoeff); 195 opticalEq->SetFactor(mass,opticalCoeff);
196 (*theEquationMap)[G4PionZero::PionZero()->G 196 (*theEquationMap)[G4PionZero::PionZero()->GetPDGEncoding()] = opticalEq;
197 197
198 opticalEq = new G4KM_OpticalEqRhs(theField, 198 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
199 mass = G4SigmaPlus::SigmaPlus()->GetPDGMass 199 mass = G4SigmaPlus::SigmaPlus()->GetPDGMass();
200 opticalCoeff = 200 opticalCoeff =
201 (*theFieldMap)[G4SigmaPlus::SigmaPlus 201 (*theFieldMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()]->GetCoeff();
202 opticalEq->SetFactor(mass,opticalCoeff); 202 opticalEq->SetFactor(mass,opticalCoeff);
203 (*theEquationMap)[G4SigmaPlus::SigmaPlus()- 203 (*theEquationMap)[G4SigmaPlus::SigmaPlus()->GetPDGEncoding()] = opticalEq;
204 204
205 opticalEq = new G4KM_OpticalEqRhs(theField, 205 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
206 mass = G4SigmaMinus::SigmaMinus()->GetPDGMa 206 mass = G4SigmaMinus::SigmaMinus()->GetPDGMass();
207 opticalCoeff = 207 opticalCoeff =
208 (*theFieldMap)[G4SigmaMinus::SigmaMin 208 (*theFieldMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()]->GetCoeff();
209 opticalEq->SetFactor(mass,opticalCoeff); 209 opticalEq->SetFactor(mass,opticalCoeff);
210 (*theEquationMap)[G4SigmaMinus::SigmaMinus( 210 (*theEquationMap)[G4SigmaMinus::SigmaMinus()->GetPDGEncoding()] = opticalEq;
211 211
212 opticalEq = new G4KM_OpticalEqRhs(theField, 212 opticalEq = new G4KM_OpticalEqRhs(theField, theNucleus);
213 mass = G4SigmaZero::SigmaZero()->GetPDGMass 213 mass = G4SigmaZero::SigmaZero()->GetPDGMass();
214 opticalCoeff = 214 opticalCoeff =
215 (*theFieldMap)[G4SigmaZero::SigmaZero 215 (*theFieldMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()]->GetCoeff();
216 opticalEq->SetFactor(mass,opticalCoeff); 216 opticalEq->SetFactor(mass,opticalCoeff);
217 (*theEquationMap)[G4SigmaZero::SigmaZero()- 217 (*theEquationMap)[G4SigmaZero::SigmaZero()->GetPDGEncoding()] = opticalEq;
218 } 218 }
219 219
220 220
221 //#define debug_1_RKPropagation 1 221 //#define debug_1_RKPropagation 1
222 //-------------------------------------------- 222 //----------------------------------------------------------------------------
223 void G4RKPropagation::Transport(G4KineticTrack 223 void G4RKPropagation::Transport(G4KineticTrackVector & active,
224 //-------------------------------------- 224 //----------------------------------------------------------------------------
225 const G4KineticTrackVector &, 225 const G4KineticTrackVector &,
226 G4double timeStep) 226 G4double timeStep)
227 { 227 {
228 // reset momentum transfer to field 228 // reset momentum transfer to field
229 theMomentumTranfer=G4ThreeVector(0,0,0); 229 theMomentumTranfer=G4ThreeVector(0,0,0);
230 230
231 // Loop over tracks 231 // Loop over tracks
232 232
233 std::vector<G4KineticTrack *>::iterator i; 233 std::vector<G4KineticTrack *>::iterator i;
234 for(i = active.begin(); i != active.end(); 234 for(i = active.begin(); i != active.end(); ++i)
235 { 235 {
236 G4double currTimeStep = timeStep; 236 G4double currTimeStep = timeStep;
237 G4KineticTrack * kt = *i; 237 G4KineticTrack * kt = *i;
238 G4int encoding = kt->GetDefinition()->Ge 238 G4int encoding = kt->GetDefinition()->GetPDGEncoding();
239 239
240 std::map <G4int, G4VNuclearField*, std:: 240 std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding);
241 241
242 G4VNuclearField* currentField=0; 242 G4VNuclearField* currentField=0;
243 if ( fieldIter != theFieldMap->end() ) c 243 if ( fieldIter != theFieldMap->end() ) currentField=fieldIter->second;
244 244
245 // debug 245 // debug
246 // if ( timeStep > 1e30 ) { 246 // if ( timeStep > 1e30 ) {
247 // G4cout << " Name :" << kt->GetDefini 247 // G4cout << " Name :" << kt->GetDefinition()->GetParticleName() << G4endl;
248 // } 248 // }
249 249
250 // Get the time of intersections with th 250 // Get the time of intersections with the nucleus surface.
251 G4double t_enter, t_leave; 251 G4double t_enter, t_leave;
252 // if the particle does not intersecate 252 // if the particle does not intersecate with the nucleus go to next particle
253 if(!GetSphereIntersectionTimes(kt, t_ent 253 if(!GetSphereIntersectionTimes(kt, t_enter, t_leave))
254 { 254 {
255 kt->SetState(G4KineticTrack::miss_nuc 255 kt->SetState(G4KineticTrack::miss_nucleus);
256 continue; 256 continue;
257 } 257 }
258 258
259 259
260 #ifdef debug_1_RKPropagation 260 #ifdef debug_1_RKPropagation
261 G4cout <<" kt,timeStep, Intersection tim 261 G4cout <<" kt,timeStep, Intersection times tenter, tleave "
262 <<kt<< " / state= " <<kt->GetState 262 <<kt<< " / state= " <<kt->GetState() <<" / " <<" "<< currTimeStep << " / " << t_enter << " / " << t_leave <<G4endl;
263 #endif 263 #endif
264 264
265 // if the particle is already outside nu 265 // if the particle is already outside nucleus go to next @@GF should never happen? check!
266 // does happen for particles added as l 266 // does happen for particles added as late....
267 // if(t_leave < 0 ) 267 // if(t_leave < 0 )
268 // { 268 // {
269 // throw G4HadronicException(__FI 269 // throw G4HadronicException(__FILE__, __LINE__, "G4RKPropagation:: Attempt to track particle past a nucleus");
270 // continue; 270 // continue;
271 // } 271 // }
272 272
273 // Apply a straight line propagation for 273 // Apply a straight line propagation for particle types
274 // not included in the model 274 // not included in the model
275 if( ! currentField ) 275 if( ! currentField )
276 { 276 {
277 if(currTimeStep == DBL_MAX)currTimeSt 277 if(currTimeStep == DBL_MAX)currTimeStep = t_leave*1.05;
278 FreeTransport(kt, currTimeStep); 278 FreeTransport(kt, currTimeStep);
279 if ( currTimeStep >= t_leave ) 279 if ( currTimeStep >= t_leave )
280 { 280 {
281 if ( kt->GetState() == G4KineticTr 281 if ( kt->GetState() == G4KineticTrack::inside )
282 { kt->SetState(G4KineticTrack::gon 282 { kt->SetState(G4KineticTrack::gone_out); }
283 else 283 else
284 { kt->SetState(G4KineticTrack::mis 284 { kt->SetState(G4KineticTrack::miss_nucleus);}
285 } else if (kt->GetState() == G4Kineti 285 } else if (kt->GetState() == G4KineticTrack::outside && currTimeStep >= t_enter ){
286 kt->SetState(G4KineticTrack::insid 286 kt->SetState(G4KineticTrack::inside);
287 } 287 }
288 288
289 continue; 289 continue;
290 } 290 }
291 291
292 if(t_enter > 0) // the particle is out. 292 if(t_enter > 0) // the particle is out. Transport free to the surface
293 { 293 {
294 if(t_enter > currTimeStep) // the pa 294 if(t_enter > currTimeStep) // the particle won't enter the nucleus
295 { 295 {
296 FreeTransport(kt, currTimeStep); 296 FreeTransport(kt, currTimeStep);
297 continue; 297 continue;
298 } 298 }
299 else 299 else
300 { 300 {
301 FreeTransport(kt, t_enter); // g 301 FreeTransport(kt, t_enter); // go to surface
302 currTimeStep -= t_enter; 302 currTimeStep -= t_enter;
303 t_leave -= t_enter; // time 303 t_leave -= t_enter; // time left to leave nucleus
304 // on the surface the particle loo 304 // on the surface the particle loose the barrier energy
305 // G4double newE = mom.e()-(*theF 305 // G4double newE = mom.e()-(*theFieldMap)[encoding]->GetBarrier();
306 // GetField = Barrier + FermiP 306 // GetField = Barrier + FermiPotential
307 G4double newE = kt->GetTrackingMom 307 G4double newE = kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition());
308 308
309 if(newE <= kt->GetActualMass()) / 309 if(newE <= kt->GetActualMass()) // the particle cannot enter the nucleus
310 { 310 {
311 // FixMe: should be "pushed bac 311 // FixMe: should be "pushed back?"
312 // for the moment take it 312 // for the moment take it past the nucleus, so we'll not worry next time..
313 FreeTransport(kt, 1.1*t_leave); 313 FreeTransport(kt, 1.1*t_leave); // take past nucleus
314 kt->SetState(G4KineticTrack::mi 314 kt->SetState(G4KineticTrack::miss_nucleus);
315 // G4cout << "G4RKPropagatio 315 // G4cout << "G4RKPropagation: Warning particle cannot enter Nucleus :" << G4endl;
316 // G4cout << " enter nucleus 316 // G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl;
317 // G4cout << " the Field "<< 317 // G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl;
318 // G4cout << " the particl 318 // G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl;
319 continue; 319 continue;
320 } 320 }
321 // 321 //
322 G4double newP = std::sqrt(newE*new 322 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
323 G4LorentzVector new4Mom(newP*kt->G 323 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
324 G4ThreeVector transfer(kt->GetTrac 324 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
325 G4ThreeVector boost= transfer / st 325 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
326 new4Mom*=G4LorentzRotation(boost); 326 new4Mom*=G4LorentzRotation(boost);
327 kt->SetTrackingMomentum(new4Mom); 327 kt->SetTrackingMomentum(new4Mom);
328 kt->SetState(G4KineticTrack::insid 328 kt->SetState(G4KineticTrack::inside);
329 329
330 /* 330 /*
331 G4cout <<" Enter Nucleus - E/Field/Sum: " 331 G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / "
332 << (*theFieldMap)[encoding]->GetField 332 << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / "
333 << kt->GetTrackingMomentum().e()-currentF 333 << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition())
334 << G4endl 334 << G4endl
335 << " Barrier / field just inside nucleus 335 << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())"
336 << (*theFieldMap)[encoding]->GetBarrier() 336 << (*theFieldMap)[encoding]->GetBarrier() << " / "
337 << (*theFieldMap)[encoding]->GetField(0.9 337 << (*theFieldMap)[encoding]->GetField(0.9999*kt->GetPosition())
338 << G4endl; 338 << G4endl;
339 */ 339 */
340 } 340 }
341 } 341 }
342 342
343 // FixMe: should I add a control on theC 343 // FixMe: should I add a control on theCutOnP here?
344 // Transport the particle into the nucle 344 // Transport the particle into the nucleus
345 // G4cerr << "RKPropagation t_leav 345 // G4cerr << "RKPropagation t_leave, curTimeStep " <<t_leave << " " <<currTimeStep<<G4endl;
346 G4bool is_exiting=false; 346 G4bool is_exiting=false;
347 if(currTimeStep > t_leave) // particle 347 if(currTimeStep > t_leave) // particle will exit from the nucleus
348 { 348 {
349 currTimeStep = t_leave; 349 currTimeStep = t_leave;
350 is_exiting=true; 350 is_exiting=true;
351 } 351 }
352 352
353 #ifdef debug_1_RKPropagation 353 #ifdef debug_1_RKPropagation
354 G4cerr << "RKPropagation is_exiting?, t_ 354 G4cerr << "RKPropagation is_exiting?, t_leave, curTimeStep " <<is_exiting<<" "<<t_leave << " " <<currTimeStep<<G4endl;
355 G4cout << "RKPropagation Ekin, field, pr 355 G4cout << "RKPropagation Ekin, field, projectile potential, p "
356 << kt->GetTrackingMomentum().e() - 356 << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
357 << kt->GetPosition()<<" " 357 << kt->GetPosition()<<" "
358 << G4endl << currentField->GetFiel 358 << G4endl << currentField->GetField(kt->GetPosition()) << " "
359 << kt->GetProjectilePotential()<< 359 << kt->GetProjectilePotential()<< G4endl
360 << kt->GetTrackingMomentum() 360 << kt->GetTrackingMomentum()
361 << G4endl; 361 << G4endl;
362 #endif 362 #endif
363 363
364 G4LorentzVector momold=kt->GetTrackingMo 364 G4LorentzVector momold=kt->GetTrackingMomentum();
365 G4ThreeVector posold=kt->GetPosition(); 365 G4ThreeVector posold=kt->GetPosition();
366 366
367 // if (currentField->GetField(kt->Get 367 // if (currentField->GetField(kt->GetPosition()) > kt->GetProjectilePotential() ||
368 if (currTimeStep > 0 && 368 if (currTimeStep > 0 &&
369 ! FieldTransport(kt, currTimeStep) 369 ! FieldTransport(kt, currTimeStep)) {
370 FreeTransport(kt,currTimeStep); 370 FreeTransport(kt,currTimeStep);
371 } 371 }
372 372
373 #ifdef debug_1_RKPropagation 373 #ifdef debug_1_RKPropagation
374 G4cout << "RKPropagation Ekin, field, p 374 G4cout << "RKPropagation Ekin, field, p "
375 << kt->GetTrackingMomentum().e() - 375 << kt->GetTrackingMomentum().e() - kt->GetTrackingMomentum().mag() << " "
376 << G4endl << currentField->GetFiel 376 << G4endl << currentField->GetField(kt->GetPosition())<< G4endl
377 << kt->GetTrackingMomentum() 377 << kt->GetTrackingMomentum()
378 << G4endl 378 << G4endl
379 << "delta p " << momold-kt->GetTra 379 << "delta p " << momold-kt->GetTrackingMomentum() << G4endl
380 << "del pos " << posold-kt->GetPos 380 << "del pos " << posold-kt->GetPosition()
381 << G4endl; 381 << G4endl;
382 #endif 382 #endif
383 383
384 // complete the transport 384 // complete the transport
385 // FixMe: in some cases there could be a 385 // FixMe: in some cases there could be a significant
386 // part to do still in the nucleu 386 // part to do still in the nucleus, or we stepped to far... depending on
387 // slope of potential 387 // slope of potential
388 G4double t_in=-1, t_out=0; // set onto 388 G4double t_in=-1, t_out=0; // set onto boundary.
389 389
390 // should go out, or are already out by 390 // should go out, or are already out by a too long step..
391 if(is_exiting || 391 if(is_exiting ||
392 (GetSphereIntersectionTimes(kt, t_ 392 (GetSphereIntersectionTimes(kt, t_in, t_out) &&t_in<0 && t_out<=0 )) // particle is exiting
393 { 393 {
394 if(t_in < 0 && t_out >= 0) //still 394 if(t_in < 0 && t_out >= 0) //still inside, transport safely out.
395 { 395 {
396 // transport free to a position th 396 // transport free to a position that is surely out of the nucleus, to avoid
397 // a new transportation and a new 397 // a new transportation and a new adding the barrier next loop.
398 G4ThreeVector savePos = kt->GetPos 398 G4ThreeVector savePos = kt->GetPosition();
399 FreeTransport(kt, t_out); 399 FreeTransport(kt, t_out);
400 // and evaluate the right the ener 400 // and evaluate the right the energy
401 G4double newE=kt->GetTrackingMomen 401 G4double newE=kt->GetTrackingMomentum().e();
402 402
403 // G4cout << " V pos/savePos << " 403 // G4cout << " V pos/savePos << "
404 // << (*theFieldMap)[encoding]- 404 // << (*theFieldMap)[encoding]->GetField(kt->GetPosition())<< " / "
405 // << (*theFieldMap)[encoding]- 405 // << (*theFieldMap)[encoding]->GetField(savePos)
406 // << G4endl; 406 // << G4endl;
407 407
408 if ( std::abs(currentField->GetFie 408 if ( std::abs(currentField->GetField(savePos)) > 0. &&
409 std::abs(currentField->GetFi 409 std::abs(currentField->GetField(kt->GetPosition())) > 0.)
410 { // FixMe GF: savePos/pos may be 410 { // FixMe GF: savePos/pos may be out of nucleus, where GetField(..)=0
411 // This wrongly adds 411 // This wrongly adds or subtracts the Barrier here while
412 // this is done later 412 // this is done later.
413 newE += currentField->GetField( 413 newE += currentField->GetField(savePos)
414 - currentField 414 - currentField->GetField(kt->GetPosition());
415 } 415 }
416 416
417 // G4cout << " go border nuc 417 // G4cout << " go border nucleus, E in/border: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
418 418
419 if(newE < kt->GetActualMass()) 419 if(newE < kt->GetActualMass())
420 { 420 {
421 #ifdef debug_1_RKPropagation 421 #ifdef debug_1_RKPropagation
422 G4cout << "RKPropagation-Transp 422 G4cout << "RKPropagation-Transport: problem with particle exiting - ignored" << G4endl;
423 G4cout << " cannot leave nucleu 423 G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
424 #endif 424 #endif
425 if (kt->GetDefinition() == G4Pr 425 if (kt->GetDefinition() == G4Proton::Proton() ||
426 kt->GetDefinition() == G4 426 kt->GetDefinition() == G4Neutron::Neutron() ) {
427 kt->SetState(G4KineticTrack: 427 kt->SetState(G4KineticTrack::captured);
428 } else { 428 } else {
429 kt->SetState(G4KineticTrack: 429 kt->SetState(G4KineticTrack::gone_out); //@@GF tofix
430 } 430 }
431 continue; // the particle canno 431 continue; // the particle cannot exit the nucleus
432 } 432 }
433 G4double newP = std::sqrt(newE*new 433 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
434 G4LorentzVector new4Mom(newP*kt->G 434 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
435 G4ThreeVector transfer(kt->GetTrac 435 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
436 G4ThreeVector boost= transfer / st 436 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
437 new4Mom*=G4LorentzRotation(boost); 437 new4Mom*=G4LorentzRotation(boost);
438 kt->SetTrackingMomentum(new4Mom); 438 kt->SetTrackingMomentum(new4Mom);
439 } 439 }
440 // add the potential barrier 440 // add the potential barrier
441 // FixMe the Coulomb field is not par 441 // FixMe the Coulomb field is not parallel to mom, this is simple approximation
442 G4double newE = kt->GetTrackingMoment 442 G4double newE = kt->GetTrackingMomentum().e()+currentField->GetField(kt->GetPosition());
443 if(newE < kt->GetActualMass()) 443 if(newE < kt->GetActualMass())
444 { // the particle cannot exit the nu 444 { // the particle cannot exit the nucleus @@@ GF check.
445 #ifdef debug_1_RKPropagation 445 #ifdef debug_1_RKPropagation
446 G4cout << " cannot leave nucleus, 446 G4cout << " cannot leave nucleus, E in/out: " << kt->GetTrackingMomentum() << " / " << newE <<G4endl;
447 #endif 447 #endif
448 if (kt->GetDefinition() == G4Proto 448 if (kt->GetDefinition() == G4Proton::Proton() ||
449 kt->GetDefinition() == G4Neu 449 kt->GetDefinition() == G4Neutron::Neutron() ) {
450 kt->SetState(G4KineticTrack::ca 450 kt->SetState(G4KineticTrack::captured);
451 } else { 451 } else {
452 kt->SetState(G4KineticTrack::go 452 kt->SetState(G4KineticTrack::gone_out); //@@GF tofix
453 } 453 }
454 continue; 454 continue;
455 } 455 }
456 G4double newP = std::sqrt(newE*newE- 456 G4double newP = std::sqrt(newE*newE- sqr(kt->GetActualMass()));
457 G4LorentzVector new4Mom(newP*kt->GetT 457 G4LorentzVector new4Mom(newP*kt->GetTrackingMomentum().vect().unit(), newE);
458 G4ThreeVector transfer(kt->GetTrackin 458 G4ThreeVector transfer(kt->GetTrackingMomentum().vect()-new4Mom.vect());
459 G4ThreeVector boost= transfer / std:: 459 G4ThreeVector boost= transfer / std::sqrt(transfer.mag2() + sqr(theNucleus->GetMass()));
460 new4Mom*=G4LorentzRotation(boost); 460 new4Mom*=G4LorentzRotation(boost);
461 kt->SetTrackingMomentum(new4Mom); 461 kt->SetTrackingMomentum(new4Mom);
462 kt->SetState(G4KineticTrack::gone_out 462 kt->SetState(G4KineticTrack::gone_out);
463 } 463 }
464 464
465 } 465 }
466 466
467 } 467 }
468 468
469 469
470 //-------------------------------------------- 470 //----------------------------------------------------------------------------
471 G4ThreeVector G4RKPropagation::GetMomentumTran 471 G4ThreeVector G4RKPropagation::GetMomentumTransfer() const
472 //-------------------------------------------- 472 //----------------------------------------------------------------------------
473 { 473 {
474 return theMomentumTranfer; 474 return theMomentumTranfer;
475 } 475 }
476 476
477 477
478 //-------------------------------------------- 478 //----------------------------------------------------------------------------
479 G4bool G4RKPropagation::FieldTransport(G4Kinet 479 G4bool G4RKPropagation::FieldTransport(G4KineticTrack * kt, const G4double timeStep)
480 //-------------------------------------------- 480 //----------------------------------------------------------------------------
481 { 481 {
482 // G4cout <<"Stepper input"<<kt->GetTrac 482 // G4cout <<"Stepper input"<<kt->GetTrackingMomentum()<<G4endl;
483 // create the integrator stepper 483 // create the integrator stepper
484 // G4Mag_EqRhs * equation = mapIter->sec 484 // G4Mag_EqRhs * equation = mapIter->second;
485 G4Mag_EqRhs * equation = (*theEquationMap)[ 485 G4Mag_EqRhs * equation = (*theEquationMap)[kt->GetDefinition()->GetPDGEncoding()];
486 G4MagIntegratorStepper * stepper = new G4Cl 486 G4MagIntegratorStepper * stepper = new G4ClassicalRK4(equation);
487 487
488 // create the integrator driver 488 // create the integrator driver
489 G4double hMin = 1.0e-25*second; // arbitr 489 G4double hMin = 1.0e-25*second; // arbitrary choice. Means 0.03 fm at c
490 G4MagInt_Driver * driver = new G4MagInt_Dri 490 G4MagInt_Driver * driver = new G4MagInt_Driver(hMin, stepper);
491 491
492 // Temporary: use driver->AccurateAdvance() 492 // Temporary: use driver->AccurateAdvance()
493 // create the G4FieldTrack needed by Accura 493 // create the G4FieldTrack needed by AccurateAdvance
494 G4double curveLength = 0; 494 G4double curveLength = 0;
495 G4FieldTrack track(kt->GetPosition(), 495 G4FieldTrack track(kt->GetPosition(),
496 kt->GetTrackingMomentum().vect().unit 496 kt->GetTrackingMomentum().vect().unit(), // momentum direction
497 curveLength, // curvelength 497 curveLength, // curvelength
498 kt->GetTrackingMomentum().e()-kt->Get 498 kt->GetTrackingMomentum().e()-kt->GetActualMass(), // kinetic energy
499 kt->GetActualMass(), // restmass 499 kt->GetActualMass(), // restmass
500 kt->GetTrackingMomentum().beta()*c_li 500 kt->GetTrackingMomentum().beta()*c_light); // velocity
501 // integrate 501 // integrate
502 G4double eps = 0.01; 502 G4double eps = 0.01;
503 // G4cout << "currTimeStep = " << currTi 503 // G4cout << "currTimeStep = " << currTimeStep << G4endl;
504 if(!driver->AccurateAdvance(track, timeStep 504 if(!driver->AccurateAdvance(track, timeStep, eps))
505 { // cannot track this particle 505 { // cannot track this particle
506 #ifdef debug_1_RKPropagation 506 #ifdef debug_1_RKPropagation
507 std::cerr << "G4RKPropagation::FieldTran 507 std::cerr << "G4RKPropagation::FieldTransport() warning: integration error."
508 << G4endl << "position " << kt->Ge 508 << G4endl << "position " << kt->GetPosition() << " 4mom " <<kt->GetTrackingMomentum()
509 <<G4endl << " timestep " <<timeSte 509 <<G4endl << " timestep " <<timeStep
510 << G4endl; 510 << G4endl;
511 #endif 511 #endif
512 delete driver; 512 delete driver;
513 delete stepper; 513 delete stepper;
514 return false; 514 return false;
515 } 515 }
516 /* 516 /*
517 G4cout <<" E/Field/Sum be4 : " <<mom.e( 517 G4cout <<" E/Field/Sum be4 : " <<mom.e() << " / "
518 << (*theFieldMap)[encoding]->GetFie 518 << (*theFieldMap)[encoding]->GetField(pos) << " / "
519 << mom.e()+(*theFieldMap)[encoding]->Ge 519 << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
520 << G4endl; 520 << G4endl;
521 */ 521 */
522 522
523 // Correct for momentum ( thus energy) tran 523 // Correct for momentum ( thus energy) transfered to nucleus, boost particle into moving nucleus frame.
524 G4ThreeVector MomentumTranfer = kt->GetTrac 524 G4ThreeVector MomentumTranfer = kt->GetTrackingMomentum().vect() - track.GetMomentum();
525 G4ThreeVector boost= MomentumTranfer / std: 525 G4ThreeVector boost= MomentumTranfer / std::sqrt (MomentumTranfer.mag2() +sqr(theNucleus->GetMass()));
526 526
527 // update the kt 527 // update the kt
528 kt->SetPosition(track.GetPosition()); 528 kt->SetPosition(track.GetPosition());
529 G4LorentzVector mom(track.GetMomentum(),std 529 G4LorentzVector mom(track.GetMomentum(),std::sqrt(track.GetMomentum().mag2() + sqr(kt->GetActualMass())));
530 mom *= G4LorentzRotation( boost ); 530 mom *= G4LorentzRotation( boost );
531 theMomentumTranfer += ( kt->GetTrackingMome 531 theMomentumTranfer += ( kt->GetTrackingMomentum() - mom ).vect();
532 kt->SetTrackingMomentum(mom); 532 kt->SetTrackingMomentum(mom);
533 533
534 // G4cout <<"Stepper output"<<kt<<" "<<k 534 // G4cout <<"Stepper output"<<kt<<" "<<kt->GetTrackingMomentum()<<" "<<kt->GetPosition()<<G4endl;
535 /* 535 /*
536 * G4ThreeVector MomentumTranfer2=kt->Get 536 * G4ThreeVector MomentumTranfer2=kt->GetTrackingMomentum().vect() - mom.vect();
537 * G4cout << " MomentumTransfer/corrected" 537 * G4cout << " MomentumTransfer/corrected" << MomentumTranfer << " " << MomentumTranfer.mag()
538 * << " " << MomentumTranfer2 << 538 * << " " << MomentumTranfer2 << " " << MomentumTranfer2.mag() << " "
539 * << MomentumTranfer-MomentumTranfer2 539 * << MomentumTranfer-MomentumTranfer2 << " "<<
540 * MomentumTranfer-MomentumTranfer2.mag 540 * MomentumTranfer-MomentumTranfer2.mag() << " " << G4endl;
541 * G4cout <<" E/Field/Sum aft : " <<mo 541 * G4cout <<" E/Field/Sum aft : " <<mom.e() << " / "
542 * << " / " << (*theFieldMap)[en 542 * << " / " << (*theFieldMap)[encoding]->GetField(pos)<< " / "
543 * << mom.e()+(*theFieldMap)[encoding] 543 * << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
544 * << G4endl; 544 * << G4endl;
545 */ 545 */
546 546
547 delete driver; 547 delete driver;
548 delete stepper; 548 delete stepper;
549 return true; 549 return true;
550 } 550 }
551 551
552 //-------------------------------------------- 552 //----------------------------------------------------------------------------
553 G4bool G4RKPropagation::FreeTransport(G4Kineti 553 G4bool G4RKPropagation::FreeTransport(G4KineticTrack * kt, const G4double timeStep)
554 //-------------------------------------------- 554 //----------------------------------------------------------------------------
555 { 555 {
556 G4ThreeVector newpos = kt->GetPosition() + 556 G4ThreeVector newpos = kt->GetPosition() +
557 timeStep*c_light/kt->GetTrackingMomen 557 timeStep*c_light/kt->GetTrackingMomentum().e() * kt->GetTrackingMomentum().vect();
558 kt->SetPosition(newpos); 558 kt->SetPosition(newpos);
559 return true; 559 return true;
560 } 560 }
561 561
562 /* 562 /*
563 G4bool G4RKPropagation::WillBeCaptured(const G 563 G4bool G4RKPropagation::WillBeCaptured(const G4KineticTrack * kt)
564 { 564 {
565 G4double radius = theOuterRadius; 565 G4double radius = theOuterRadius;
566 566
567 // evaluate the final energy. Il will be captu 567 // evaluate the final energy. Il will be captured if newE or newP < 0
568 G4ParticleDefinition * definition = kt->GetD 568 G4ParticleDefinition * definition = kt->GetDefinition();
569 G4double mass = definition->GetPDGMass(); 569 G4double mass = definition->GetPDGMass();
570 G4ThreeVector pos = kt->GetPosition(); 570 G4ThreeVector pos = kt->GetPosition();
571 G4LorentzVector mom = kt->GetTrackingMomentu 571 G4LorentzVector mom = kt->GetTrackingMomentum();
572 G4VNuclearField * field = (*theFieldMap)[def 572 G4VNuclearField * field = (*theFieldMap)[definition->GetPDGEncoding()];
573 G4ThreeVector newPos(0, 0, radius); // to ge 573 G4ThreeVector newPos(0, 0, radius); // to get the field on the surface
574 574
575 G4double newE = mom.e()+field->GetField(pos) 575 G4double newE = mom.e()+field->GetField(pos)-field->GetField(newPos);
576 576
577 return ((newE < mass) ? false : true); 577 return ((newE < mass) ? false : true);
578 } 578 }
579 */ 579 */
580 580
581 581
582 582
583 //-------------------------------------------- 583 //----------------------------------------------------------------------------
584 G4bool G4RKPropagation::GetSphereIntersectionT 584 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4double radius,
585 //-------------------------------------- 585 //----------------------------------------------------------------------------
586 const G4ThreeVector & currentPos, 586 const G4ThreeVector & currentPos,
587 const G4LorentzVector & momentum, 587 const G4LorentzVector & momentum,
588 G4double & t1, G4double & t2) 588 G4double & t1, G4double & t2)
589 { 589 {
590 G4ThreeVector speed = momentum.vect()/momen 590 G4ThreeVector speed = momentum.vect()/momentum.e(); // boost vector
591 G4double scalarProd = currentPos.dot(speed) 591 G4double scalarProd = currentPos.dot(speed);
592 G4double speedMag2 = speed.mag2(); 592 G4double speedMag2 = speed.mag2();
593 G4double sqrtArg = scalarProd*scalarProd - 593 G4double sqrtArg = scalarProd*scalarProd -
594 speedMag2*(currentPos.mag2()-radius*r 594 speedMag2*(currentPos.mag2()-radius*radius);
595 if(sqrtArg <= 0.) // particle will not inte 595 if(sqrtArg <= 0.) // particle will not intersect the sphere
596 { 596 {
597 // G4cout << " GetSphereIntersection 597 // G4cout << " GetSphereIntersectionTimes sqrtArg negative: " << sqrtArg << G4endl;
598 return false; 598 return false;
599 } 599 }
600 t1 = (-scalarProd - std::sqrt(sqrtArg))/spe 600 t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
601 t2 = (-scalarProd + std::sqrt(sqrtArg))/spe 601 t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
602 return true; 602 return true;
603 } 603 }
604 604
605 //-------------------------------------------- 605 //----------------------------------------------------------------------------
606 G4bool G4RKPropagation::GetSphereIntersectionT 606 G4bool G4RKPropagation::GetSphereIntersectionTimes(const G4KineticTrack * kt,
607 G4double & t1, G4double & t2) 607 G4double & t1, G4double & t2)
608 { 608 {
609 G4double radius = theOuterRadius + 3*fermi; 609 G4double radius = theOuterRadius + 3*fermi; // "safety" of 3 fermi
610 G4ThreeVector speed = kt->GetTrackingMoment 610 G4ThreeVector speed = kt->GetTrackingMomentum().vect()/kt->GetTrackingMomentum().e(); // bost vector
611 G4double scalarProd = kt->GetPosition().dot 611 G4double scalarProd = kt->GetPosition().dot(speed);
612 G4double speedMag2 = speed.mag2(); 612 G4double speedMag2 = speed.mag2();
613 G4double sqrtArg = scalarProd*scalarProd - 613 G4double sqrtArg = scalarProd*scalarProd -
614 speedMag2*(kt->GetPosition().mag2()-r 614 speedMag2*(kt->GetPosition().mag2()-radius*radius);
615 if(sqrtArg <= 0.) // particle will not inte 615 if(sqrtArg <= 0.) // particle will not intersect the sphere
616 { 616 {
617 return false; 617 return false;
618 } 618 }
619 t1 = (-scalarProd - std::sqrt(sqrtArg))/spe 619 t1 = (-scalarProd - std::sqrt(sqrtArg))/speedMag2/c_light;
620 t2 = (-scalarProd + std::sqrt(sqrtArg))/spe 620 t2 = (-scalarProd + std::sqrt(sqrtArg))/speedMag2/c_light;
621 return true; 621 return true;
622 } 622 }
623 623
624 // Implementation methods 624 // Implementation methods
625 625
626 //-------------------------------------------- 626 //----------------------------------------------------------------------------
627 void G4RKPropagation::delete_FieldsAndMap( 627 void G4RKPropagation::delete_FieldsAndMap(
628 //-------------------------------------- 628 //----------------------------------------------------------------------------
629 std::map <G4int, G4VNuclearField *, std: 629 std::map <G4int, G4VNuclearField *, std::less<G4int> > * aMap)
630 { 630 {
631 if(aMap) 631 if(aMap)
632 { 632 {
633 std::map <G4int, G4VNuclearField *, std: 633 std::map <G4int, G4VNuclearField *, std::less<G4int> >::iterator cur;
634 for(cur = aMap->begin(); cur != aMap->en 634 for(cur = aMap->begin(); cur != aMap->end(); ++cur)
635 delete (*cur).second; 635 delete (*cur).second;
636 636
637 aMap->clear(); 637 aMap->clear();
638 delete aMap; 638 delete aMap;
639 } 639 }
640 640
641 } 641 }
642 642
643 //-------------------------------------------- 643 //----------------------------------------------------------------------------
644 void G4RKPropagation::delete_EquationsAndMap( 644 void G4RKPropagation::delete_EquationsAndMap(
645 //-------------------------------------- 645 //----------------------------------------------------------------------------
646 std::map <G4int, G4Mag_EqRhs *, std::les 646 std::map <G4int, G4Mag_EqRhs *, std::less<G4int> > * aMap)
647 { 647 {
648 if(aMap) 648 if(aMap)
649 { 649 {
650 std::map <G4int, G4Mag_EqRhs *, std::les 650 std::map <G4int, G4Mag_EqRhs *, std::less<G4int> >::iterator cur;
651 for(cur = aMap->begin(); cur != aMap->en 651 for(cur = aMap->begin(); cur != aMap->end(); ++cur)
652 delete (*cur).second; 652 delete (*cur).second;
653 653
654 aMap->clear(); 654 aMap->clear();
655 delete aMap; 655 delete aMap;
656 } 656 }
657 } 657 }
658 658