| 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 // $Id$
26 // 27 //
27 /// \brief This class is a slightly modified v 28 /// \brief This class is a slightly modified version of G4Transportation
28 /// initially written by John Apostolak 29 /// initially written by John Apostolakis and colleagues
29 /// But it should use the exact same al 30 /// But it should use the exact same algorithm
30 // 31 //
31 // Contact : Mathieu Karamitros (kara (AT) cen 32 // Contact : Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
32 // 33 //
33 // History : 34 // History :
34 // ----------- 35 // -----------
35 // =========================================== 36 // =======================================================================
36 // Modified: 37 // Modified:
37 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravi 38 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment
38 // 20 Nov 2008, J.Apostolakis: Push safety 39 // 20 Nov 2008, J.Apostolakis: Push safety to helper - after ComputeSafety
39 // 9 Nov 2007, J.Apostolakis: Flag for sho 40 // 9 Nov 2007, J.Apostolakis: Flag for short steps, push safety to helper
40 // 19 Jan 2006, P.MoraDeFreitas: Fix for su 41 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
41 // 11 Aug 2004, M.Asai: Add G4VSensitiveDet 42 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
42 // 21 June 2003, J.Apostolakis: Calling fiel 43 // 21 June 2003, J.Apostolakis: Calling field manager with
43 // track, to enable it to 44 // track, to enable it to configure its accuracy
44 // 13 May 2003, J.Apostolakis: Zero field a 45 // 13 May 2003, J.Apostolakis: Zero field areas now taken into
45 // account correclty in al 46 // account correclty in all cases (thanks to W Pokorski).
46 // 29 June 2001, J.Apostolakis, D.Cote-Ahern 47 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
47 // correction for spin tra 48 // correction for spin tracking
48 // 20 Febr 2001, J.Apostolakis: update for 49 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack
49 // 22 Sept 2000, V.Grichine: update of K 50 // 22 Sept 2000, V.Grichine: update of Kinetic Energy
50 // ------------------------------------------ 51 // ---------------------------------------------------
51 // 10 Oct 2011, M.Karamitros: G4ITTranspo 52 // 10 Oct 2011, M.Karamitros: G4ITTransportation created
52 // Created: 19 March 1997, J. Apostolakis 53 // Created: 19 March 1997, J. Apostolakis
53 // =========================================== 54 // =======================================================================
54 // 55 //
55 // ------------------------------------------- 56 // -------------------------------------------------------------------
56 57
57 #include "G4ITTransportation.hh" 58 #include "G4ITTransportation.hh"
58 <<
59 #include "G4ChordFinder.hh" <<
60 #include "G4FieldManager.hh" <<
61 #include "G4FieldManagerStore.hh" <<
62 #include "G4IT.hh" 59 #include "G4IT.hh"
63 #include "G4ITNavigator.hh" << 60 #include "G4TrackingInformation.hh"
64 #include "G4ITSafetyHelper.hh" << 61 #include "G4SystemOfUnits.hh"
>> 62 #include "G4TransportationManager.hh"
65 #include "G4ITTransportationManager.hh" 63 #include "G4ITTransportationManager.hh"
66 #include "G4LowEnergyEmProcessSubType.hh" <<
67 #include "G4ParticleTable.hh" <<
68 #include "G4ProductionCutsTable.hh" 64 #include "G4ProductionCutsTable.hh"
>> 65 #include "G4ParticleTable.hh"
>> 66 #include "G4ITNavigator.hh"
69 #include "G4PropagatorInField.hh" 67 #include "G4PropagatorInField.hh"
70 #include "G4ReferenceCast.hh" << 68 #include "G4FieldManager.hh"
71 #include "G4SystemOfUnits.hh" << 69 #include "G4ChordFinder.hh"
72 #include "G4TrackingInformation.hh" << 70 #include "G4ITSafetyHelper.hh"
73 #include "G4TransportationManager.hh" << 71 #include "G4FieldManagerStore.hh"
74 #include "G4UnitsTable.hh" <<
75 72
76 #include <memory> << 73 #include "G4UnitsTable.hh"
>> 74 #include "G4ReferenceCast.hh"
77 75
78 class G4VSensitiveDetector; 76 class G4VSensitiveDetector;
79 77
80 #ifndef PrepareState 78 #ifndef PrepareState
81 # define PrepareState() << 79 #define PrepareState() G4ITTransportationState* __state = this->GetState<G4ITTransportationState>();
82 G4ITTransportationState* __state = this->G <<
83 #endif 80 #endif
84 81
85 #ifndef State 82 #ifndef State
86 #define State(theXInfo) (__state->theXInfo) 83 #define State(theXInfo) (__state->theXInfo)
87 #endif 84 #endif
88 85
89 //#define DEBUG_MEM 86 //#define DEBUG_MEM
90 87
91 #ifdef DEBUG_MEM 88 #ifdef DEBUG_MEM
92 #include "G4MemStat.hh" 89 #include "G4MemStat.hh"
93 using namespace G4MemStat; 90 using namespace G4MemStat;
94 using G4MemStat::MemStat; 91 using G4MemStat::MemStat;
95 #endif 92 #endif
96 93
97 //#define G4DEBUG_TRANSPORT 1 94 //#define G4DEBUG_TRANSPORT 1
98 95
99 G4ITTransportation::G4ITTransportation(const G 96 G4ITTransportation::G4ITTransportation(const G4String& aName, int verbose) :
100 G4VITProcess(aName, fTransportation), 97 G4VITProcess(aName, fTransportation),
101 fThreshold_Warning_Energy(100 * MeV), 98 fThreshold_Warning_Energy(100 * MeV),
102 fThreshold_Important_Energy(250 * MeV), 99 fThreshold_Important_Energy(250 * MeV),
103 << 100 fThresholdTrials(10),
104 fUnimportant_Energy(1 * MeV), // Old defa << 101 fUnimportant_Energy(1 * MeV), // Not used
>> 102 fSumEnergyKilled(0.0),
>> 103 fMaxEnergyKilled(0.0),
>> 104 fShortStepOptimisation(false), // Old default: true (=fast short steps)
105 fVerboseLevel(verbose) 105 fVerboseLevel(verbose)
106 { 106 {
107 pParticleChange = &fParticleChange; 107 pParticleChange = &fParticleChange;
108 G4TransportationManager* transportMgr; 108 G4TransportationManager* transportMgr;
109 transportMgr = G4TransportationManager::GetT 109 transportMgr = G4TransportationManager::GetTransportationManager();
110 G4ITTransportationManager* ITtransportMgr; 110 G4ITTransportationManager* ITtransportMgr;
111 ITtransportMgr = G4ITTransportationManager:: 111 ITtransportMgr = G4ITTransportationManager::GetTransportationManager();
112 fLinearNavigator = ITtransportMgr->GetNaviga 112 fLinearNavigator = ITtransportMgr->GetNavigatorForTracking();
113 fFieldPropagator = transportMgr->GetPropagat 113 fFieldPropagator = transportMgr->GetPropagatorInField();
114 fpSafetyHelper = ITtransportMgr->GetSafetyHe 114 fpSafetyHelper = ITtransportMgr->GetSafetyHelper(); // New
115 115
116 // Cannot determine whether a field exists h 116 // Cannot determine whether a field exists here, as it would
117 // depend on the relative order of creating 117 // depend on the relative order of creating the detector's
118 // field and this process. That order is no 118 // field and this process. That order is not guaranted.
119 // Instead later the method DoesGlobalFieldE 119 // Instead later the method DoesGlobalFieldExist() is called
120 120
121 enableAtRestDoIt = false; 121 enableAtRestDoIt = false;
122 enableAlongStepDoIt = true; 122 enableAlongStepDoIt = true;
123 enablePostStepDoIt = true; 123 enablePostStepDoIt = true;
124 SetProcessSubType(fLowEnergyTransportation); << 124 SetProcessSubType(60);
125 SetInstantiateProcessState(true); 125 SetInstantiateProcessState(true);
126 G4VITProcess::SetInstantiateProcessState(fal 126 G4VITProcess::SetInstantiateProcessState(false);
127 fInstantiateProcessState = true; 127 fInstantiateProcessState = true;
128 128
129 G4VITProcess::fpState = std::make_shared<G4I << 129 G4VITProcess::fpState.reset(new G4ITTransportationState());
130 /* 130 /*
131 if(fTransportationState == 0) 131 if(fTransportationState == 0)
132 { 132 {
133 G4cout << "KILL in G4ITTransportation::G4IT 133 G4cout << "KILL in G4ITTransportation::G4ITTransportation" << G4endl;
134 abort(); 134 abort();
135 } 135 }
136 */ 136 */
137 } 137 }
138 138
139 G4ITTransportation::G4ITTransportation(const G 139 G4ITTransportation::G4ITTransportation(const G4ITTransportation& right) :
140 G4VITProcess(right) 140 G4VITProcess(right)
141 { 141 {
142 // Copy attributes 142 // Copy attributes
143 fVerboseLevel = right.fVerboseLevel; 143 fVerboseLevel = right.fVerboseLevel;
144 fThreshold_Warning_Energy = right.fThreshold 144 fThreshold_Warning_Energy = right.fThreshold_Warning_Energy;
145 fThreshold_Important_Energy = right.fThresho 145 fThreshold_Important_Energy = right.fThreshold_Important_Energy;
146 fThresholdTrials = right.fThresholdTrials; 146 fThresholdTrials = right.fThresholdTrials;
147 fUnimportant_Energy = right.fUnimportant_Ene 147 fUnimportant_Energy = right.fUnimportant_Energy;
148 fSumEnergyKilled = right.fSumEnergyKilled; 148 fSumEnergyKilled = right.fSumEnergyKilled;
149 fMaxEnergyKilled = right.fMaxEnergyKilled; 149 fMaxEnergyKilled = right.fMaxEnergyKilled;
150 fShortStepOptimisation = right.fShortStepOpt 150 fShortStepOptimisation = right.fShortStepOptimisation;
151 151
152 // Setup Navigators 152 // Setup Navigators
153 G4TransportationManager* transportMgr; 153 G4TransportationManager* transportMgr;
154 transportMgr = G4TransportationManager::GetT 154 transportMgr = G4TransportationManager::GetTransportationManager();
155 G4ITTransportationManager* ITtransportMgr; 155 G4ITTransportationManager* ITtransportMgr;
156 ITtransportMgr = G4ITTransportationManager:: 156 ITtransportMgr = G4ITTransportationManager::GetTransportationManager();
157 fLinearNavigator = ITtransportMgr->GetNaviga 157 fLinearNavigator = ITtransportMgr->GetNavigatorForTracking();
158 fFieldPropagator = transportMgr->GetPropagat 158 fFieldPropagator = transportMgr->GetPropagatorInField();
159 fpSafetyHelper = ITtransportMgr->GetSafetyHe 159 fpSafetyHelper = ITtransportMgr->GetSafetyHelper(); // New
160 160
161 // Cannot determine whether a field exists h 161 // Cannot determine whether a field exists here, as it would
162 // depend on the relative order of creating 162 // depend on the relative order of creating the detector's
163 // field and this process. That order is no 163 // field and this process. That order is not guaranted.
164 // Instead later the method DoesGlobalFieldE 164 // Instead later the method DoesGlobalFieldExist() is called
165 165
166 enableAtRestDoIt = false; 166 enableAtRestDoIt = false;
167 enableAlongStepDoIt = true; 167 enableAlongStepDoIt = true;
168 enablePostStepDoIt = true; 168 enablePostStepDoIt = true;
169 169
170 pParticleChange = &fParticleChange; 170 pParticleChange = &fParticleChange;
171 SetInstantiateProcessState(true); 171 SetInstantiateProcessState(true);
172 G4VITProcess::SetInstantiateProcessState(fal 172 G4VITProcess::SetInstantiateProcessState(false);
173 fInstantiateProcessState = right.fInstantiat 173 fInstantiateProcessState = right.fInstantiateProcessState;
174 } 174 }
175 175
176 G4ITTransportation& G4ITTransportation::operat 176 G4ITTransportation& G4ITTransportation::operator=(const G4ITTransportation& /*right*/)
177 { 177 {
178 // if (this == &right) return *this; 178 // if (this == &right) return *this;
179 return *this; 179 return *this;
180 } 180 }
181 181
182 ////////////////////////////////////////////// 182 //////////////////////////////////////////////////////////////////////////////
183 /// Process State 183 /// Process State
184 ////////////////////////////////////////////// 184 //////////////////////////////////////////////////////////////////////////////
185 G4ITTransportation::G4ITTransportationState::G 185 G4ITTransportation::G4ITTransportationState::G4ITTransportationState() :
186 fCurrentTouchableHandle(nullptr) << 186 G4ProcessState(), fCurrentTouchableHandle(0)
187 { 187 {
188 fTransportEndPosition = G4ThreeVector(0, 0, 188 fTransportEndPosition = G4ThreeVector(0, 0, 0);
189 fTransportEndMomentumDir = G4ThreeVector(0, 189 fTransportEndMomentumDir = G4ThreeVector(0, 0, 0);
190 fTransportEndKineticEnergy = -1; 190 fTransportEndKineticEnergy = -1;
191 fTransportEndSpin = G4ThreeVector(0, 0, 0); 191 fTransportEndSpin = G4ThreeVector(0, 0, 0);
192 fMomentumChanged = false; 192 fMomentumChanged = false;
193 fEnergyChanged = false; 193 fEnergyChanged = false;
194 fEndGlobalTimeComputed = false; 194 fEndGlobalTimeComputed = false;
195 fCandidateEndGlobalTime = -1; 195 fCandidateEndGlobalTime = -1;
196 fParticleIsLooping = false; 196 fParticleIsLooping = false;
197 197
198 static G4ThreadLocal G4TouchableHandle *null << 198 static G4ThreadLocal G4TouchableHandle *nullTouchableHandle = 0;
199 if (nullTouchableHandle == nullptr) nullTouc << 199 if (!nullTouchableHandle) nullTouchableHandle = new G4TouchableHandle;
200 // Points to (G4VTouchable*) 0 200 // Points to (G4VTouchable*) 0
201 201
202 fCurrentTouchableHandle = *nullTouchableHand 202 fCurrentTouchableHandle = *nullTouchableHandle;
203 fGeometryLimitedStep = false; 203 fGeometryLimitedStep = false;
204 fPreviousSftOrigin = G4ThreeVector(0, 0, 0); 204 fPreviousSftOrigin = G4ThreeVector(0, 0, 0);
205 fPreviousSafety = 0.0; 205 fPreviousSafety = 0.0;
206 fNoLooperTrials = 0; << 206 fNoLooperTrials = false;
207 fEndPointDistance = -1; 207 fEndPointDistance = -1;
208 } 208 }
209 209
210 G4ITTransportation::G4ITTransportationState::~ 210 G4ITTransportation::G4ITTransportationState::~G4ITTransportationState()
211 { 211 {
212 ; 212 ;
213 } 213 }
214 214
215 G4ITTransportation::~G4ITTransportation() 215 G4ITTransportation::~G4ITTransportation()
216 { 216 {
217 #ifdef G4VERBOSE 217 #ifdef G4VERBOSE
218 if ((fVerboseLevel > 0) && (fSumEnergyKilled 218 if ((fVerboseLevel > 0) && (fSumEnergyKilled > 0.0))
219 { 219 {
220 G4cout << " G4ITTransportation: Statistics 220 G4cout << " G4ITTransportation: Statistics for looping particles "
221 << G4endl; 221 << G4endl;
222 G4cout << " Sum of energy of loopers kil 222 G4cout << " Sum of energy of loopers killed: "
223 << fSumEnergyKilled << G4endl; 223 << fSumEnergyKilled << G4endl;
224 G4cout << " Max energy of loopers killed 224 G4cout << " Max energy of loopers killed: "
225 << fMaxEnergyKilled << G4endl; 225 << fMaxEnergyKilled << G4endl;
226 } 226 }
227 #endif 227 #endif
228 } 228 }
229 229
230 void G4ITTransportation::BuildPhysicsTable(con 230 void G4ITTransportation::BuildPhysicsTable(const G4ParticleDefinition&)
231 { 231 {
232 fpSafetyHelper->InitialiseHelper(); 232 fpSafetyHelper->InitialiseHelper();
233 } 233 }
234 234
235 G4bool G4ITTransportation::DoesGlobalFieldExis 235 G4bool G4ITTransportation::DoesGlobalFieldExist()
236 { 236 {
237 G4TransportationManager* transportMgr; 237 G4TransportationManager* transportMgr;
238 transportMgr = G4TransportationManager::GetT 238 transportMgr = G4TransportationManager::GetTransportationManager();
239 239
240 // fFieldExists= transportMgr->GetFieldManag 240 // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist();
241 // return fFieldExists; 241 // return fFieldExists;
242 return transportMgr->GetFieldManager()->Does 242 return transportMgr->GetFieldManager()->DoesFieldExist();
243 } 243 }
244 244
245 ////////////////////////////////////////////// 245 //////////////////////////////////////////////////////////////////////////
246 // 246 //
247 // Responsibilities: 247 // Responsibilities:
248 // Find whether the geometry limits the Ste 248 // Find whether the geometry limits the Step, and to what length
249 // Calculate the new value of the safety an 249 // Calculate the new value of the safety and return it.
250 // Store the final time, position and momen 250 // Store the final time, position and momentum.
251 G4double 251 G4double
252 G4ITTransportation:: 252 G4ITTransportation::
253 AlongStepGetPhysicalInteractionLength(const G4 253 AlongStepGetPhysicalInteractionLength(const G4Track& track,
254 G4double 254 G4double,
255 G4double 255 G4double currentMinimumStep,
256 G4double 256 G4double& currentSafety,
257 G4GPILSe 257 G4GPILSelection* selection)
258 { 258 {
259 PrepareState(); << 259 PrepareState()
260 G4double geometryStepLength(-1.0), newSafety 260 G4double geometryStepLength(-1.0), newSafety(-1.0);
261 261
262 State(fParticleIsLooping) = false; 262 State(fParticleIsLooping) = false;
263 State(fEndGlobalTimeComputed) = false; 263 State(fEndGlobalTimeComputed) = false;
264 State(fGeometryLimitedStep) = false; 264 State(fGeometryLimitedStep) = false;
265 265
266 // Initial actions moved to StartTrack() 266 // Initial actions moved to StartTrack()
267 // -------------------------------------- 267 // --------------------------------------
268 // Note: in case another process changes tou 268 // Note: in case another process changes touchable handle
269 // it will be necessary to add here (for 269 // it will be necessary to add here (for all steps)
270 // State(fCurrentTouchableHandle) = track.Ge 270 // State(fCurrentTouchableHandle) = track.GetTouchableHandle();
271 271
272 // GPILSelection is set to defaule value of 272 // GPILSelection is set to defaule value of CandidateForSelection
273 // It is a return value 273 // It is a return value
274 // 274 //
275 *selection = CandidateForSelection; 275 *selection = CandidateForSelection;
276 276
277 // Get initial Energy/Momentum of the track 277 // Get initial Energy/Momentum of the track
278 // 278 //
279 const G4DynamicParticle* pParticle = track.G 279 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
280 // const G4ParticleDefinition* pParticleDef 280 // const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
281 G4ThreeVector startMomentumDir = pParticle-> 281 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection();
282 G4ThreeVector startPosition = track.GetPosit 282 G4ThreeVector startPosition = track.GetPosition();
283 283
284 // G4double theTime = track.GetGlob 284 // G4double theTime = track.GetGlobalTime() ;
285 285
286 // The Step Point safety can be limited by o 286 // The Step Point safety can be limited by other geometries and/or the
287 // assumptions of any process - it's not alw 287 // assumptions of any process - it's not always the geometrical safety.
288 // We calculate the starting point's isotrop 288 // We calculate the starting point's isotropic safety here.
289 // 289 //
290 G4ThreeVector OriginShift = startPosition - 290 G4ThreeVector OriginShift = startPosition - State(fPreviousSftOrigin);
291 G4double MagSqShift = OriginShift.mag2(); 291 G4double MagSqShift = OriginShift.mag2();
292 if (MagSqShift >= sqr(State(fPreviousSafety) 292 if (MagSqShift >= sqr(State(fPreviousSafety)))
293 { 293 {
294 currentSafety = 0.0; 294 currentSafety = 0.0;
295 } 295 }
296 else 296 else
297 { 297 {
298 currentSafety = State(fPreviousSafety) - s 298 currentSafety = State(fPreviousSafety) - std::sqrt(MagSqShift);
299 } 299 }
300 300
301 // Is the particle charged ? 301 // Is the particle charged ?
302 // 302 //
303 G4double particleCharge = pParticle->GetChar 303 G4double particleCharge = pParticle->GetCharge();
304 304
305 // There is no need to locate the current vo 305 // There is no need to locate the current volume. It is Done elsewhere:
306 // On track construction 306 // On track construction
307 // By the tracking, after all AlongStepDoI 307 // By the tracking, after all AlongStepDoIts, in "Relocation"
308 308
309 // Check whether the particle have an (EM) f 309 // Check whether the particle have an (EM) field force exerting upon it
310 // 310 //
311 G4FieldManager* fieldMgr = nullptr; << 311 G4FieldManager* fieldMgr = 0;
312 G4bool fieldExertsForce = false; 312 G4bool fieldExertsForce = false;
313 if ((particleCharge != 0.0)) 313 if ((particleCharge != 0.0))
314 { 314 {
315 fieldMgr = fFieldPropagator->FindAndSetFie 315 fieldMgr = fFieldPropagator->FindAndSetFieldManager(track.GetVolume());
316 if (fieldMgr != nullptr) << 316 if (fieldMgr != 0)
317 { 317 {
318 // Message the field Manager, to configu 318 // Message the field Manager, to configure it for this track
319 fieldMgr->ConfigureForTrack(&track); 319 fieldMgr->ConfigureForTrack(&track);
320 // Moved here, in order to allow a trans 320 // Moved here, in order to allow a transition
321 // from a zero-field status (with fie 321 // from a zero-field status (with fieldMgr->(field)0
322 // to a finite field status 322 // to a finite field status
323 323
324 // If the field manager has no field, th 324 // If the field manager has no field, there is no field !
325 fieldExertsForce = (fieldMgr->GetDetecto << 325 fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
326 } 326 }
327 } 327 }
328 328
329 // G4cout << " G4Transport: field exerts fo 329 // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
330 // << " fieldMgr= " << fieldMgr << G4endl 330 // << " fieldMgr= " << fieldMgr << G4endl;
331 331
332 // Choose the calculation of the transportat 332 // Choose the calculation of the transportation: Field or not
333 // 333 //
334 if (!fieldExertsForce) 334 if (!fieldExertsForce)
335 { 335 {
336 G4double linearStepLength; 336 G4double linearStepLength;
337 if (fShortStepOptimisation && (currentMini 337 if (fShortStepOptimisation && (currentMinimumStep <= currentSafety))
338 { 338 {
339 // The Step is guaranteed to be taken 339 // The Step is guaranteed to be taken
340 // 340 //
341 geometryStepLength = currentMinimumStep; 341 geometryStepLength = currentMinimumStep;
342 State(fGeometryLimitedStep) = false; 342 State(fGeometryLimitedStep) = false;
343 } 343 }
344 else 344 else
345 { 345 {
346 // Find whether the straight path inter 346 // Find whether the straight path intersects a volume
347 // 347 //
348 // fLinearNavigator->SetNavigatorState(G 348 // fLinearNavigator->SetNavigatorState(GetIT(track)->GetTrackingInfo()->GetNavigatorState());
349 linearStepLength = fLinearNavigator->Com 349 linearStepLength = fLinearNavigator->ComputeStep(startPosition,
350 350 startMomentumDir,
351 351 currentMinimumStep,
352 352 newSafety);
353 353
354 // G4cout << "linearStepLength : " << G4 354 // G4cout << "linearStepLength : " << G4BestUnit(linearStepLength,"Length")
355 // << " | currentMinimumStep: " << cu 355 // << " | currentMinimumStep: " << currentMinimumStep
356 // << " | trackID: " << track.GetTrac 356 // << " | trackID: " << track.GetTrackID() << G4endl;
357 357
358 // Remember last safety origin & value. 358 // Remember last safety origin & value.
359 // 359 //
360 State(fPreviousSftOrigin) = startPositio 360 State(fPreviousSftOrigin) = startPosition;
361 State(fPreviousSafety) = newSafety; 361 State(fPreviousSafety) = newSafety;
362 362
363 G4TrackStateManager& trackStateMan = Get 363 G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo()
364 ->GetTrackStateManager(); 364 ->GetTrackStateManager();
365 fpSafetyHelper->LoadTrackState(trackStat 365 fpSafetyHelper->LoadTrackState(trackStateMan);
366 // fpSafetyHelper->SetTrackState(state); 366 // fpSafetyHelper->SetTrackState(state);
367 fpSafetyHelper->SetCurrentSafety(newSafe 367 fpSafetyHelper->SetCurrentSafety(newSafety,
368 State(f 368 State(fTransportEndPosition));
369 fpSafetyHelper->ResetTrackState(); 369 fpSafetyHelper->ResetTrackState();
370 370
371 // The safety at the initial point has b 371 // The safety at the initial point has been re-calculated:
372 // 372 //
373 currentSafety = newSafety; 373 currentSafety = newSafety;
374 374
375 State(fGeometryLimitedStep) = (linearSte 375 State(fGeometryLimitedStep) = (linearStepLength <= currentMinimumStep);
376 if (State(fGeometryLimitedStep)) 376 if (State(fGeometryLimitedStep))
377 { 377 {
378 // The geometry limits the Step size ( 378 // The geometry limits the Step size (an intersection was found.)
379 geometryStepLength = linearStepLength; 379 geometryStepLength = linearStepLength;
380 } 380 }
381 else 381 else
382 { 382 {
383 // The full Step is taken. 383 // The full Step is taken.
384 geometryStepLength = currentMinimumSte 384 geometryStepLength = currentMinimumStep;
385 } 385 }
386 } 386 }
387 State(fEndPointDistance) = geometryStepLen 387 State(fEndPointDistance) = geometryStepLength;
388 388
389 // Calculate final position 389 // Calculate final position
390 // 390 //
391 State(fTransportEndPosition) = startPositi 391 State(fTransportEndPosition) = startPosition
392 + geometryStepLength * startMomentumDi 392 + geometryStepLength * startMomentumDir;
393 393
394 // Momentum direction, energy and polarisa 394 // Momentum direction, energy and polarisation are unchanged by transport
395 // 395 //
396 State(fTransportEndMomentumDir) = startMom 396 State(fTransportEndMomentumDir) = startMomentumDir;
397 State(fTransportEndKineticEnergy) = track. 397 State(fTransportEndKineticEnergy) = track.GetKineticEnergy();
398 State(fTransportEndSpin) = track.GetPolari 398 State(fTransportEndSpin) = track.GetPolarization();
399 State(fParticleIsLooping) = false; 399 State(fParticleIsLooping) = false;
400 State(fMomentumChanged) = false; 400 State(fMomentumChanged) = false;
401 State(fEndGlobalTimeComputed) = true; 401 State(fEndGlobalTimeComputed) = true;
402 State(theInteractionTimeLeft) = State(fEnd 402 State(theInteractionTimeLeft) = State(fEndPointDistance)
403 / track.GetVelocity(); 403 / track.GetVelocity();
404 State(fCandidateEndGlobalTime) = State(the 404 State(fCandidateEndGlobalTime) = State(theInteractionTimeLeft)
405 + track.GetGlobalTime(); 405 + track.GetGlobalTime();
406 /* 406 /*
407 G4cout << "track.GetVelocity() : " 407 G4cout << "track.GetVelocity() : "
408 << track.GetVelocity() << G4endl; 408 << track.GetVelocity() << G4endl;
409 G4cout << "State(endpointDistance) : " 409 G4cout << "State(endpointDistance) : "
410 << G4BestUnit(State(endpointDistanc 410 << G4BestUnit(State(endpointDistance),"Length") << G4endl;
411 G4cout << "State(theInteractionTimeLeft) : 411 G4cout << "State(theInteractionTimeLeft) : "
412 << G4BestUnit(State(theInteractionT 412 << G4BestUnit(State(theInteractionTimeLeft),"Time") << G4endl;
413 G4cout << "track.GetGlobalTime() : " 413 G4cout << "track.GetGlobalTime() : "
414 << G4BestUnit(track.GetGlobalTime() 414 << G4BestUnit(track.GetGlobalTime(),"Time") << G4endl;
415 */ 415 */
416 } 416 }
417 else // A field exerts force 417 else // A field exerts force
418 { 418 {
419 419
420 G4ExceptionDescription exceptionDescriptio 420 G4ExceptionDescription exceptionDescription;
421 exceptionDescription 421 exceptionDescription
422 << "ITTransportation does not support 422 << "ITTransportation does not support external fields.";
423 exceptionDescription 423 exceptionDescription
424 << " If you are dealing with a tradiat 424 << " If you are dealing with a tradiational MC simulation, ";
425 exceptionDescription << "please use G4Tran 425 exceptionDescription << "please use G4Transportation.";
426 426
427 G4Exception("G4ITTransportation::AlongStep 427 G4Exception("G4ITTransportation::AlongStepGetPhysicalInteractionLength",
428 "NoExternalFieldSupport", Fata 428 "NoExternalFieldSupport", FatalException, exceptionDescription);
429 /* 429 /*
430 G4double momentumMagnitude = pParti 430 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
431 // G4ThreeVector EndUnitMomentum 431 // G4ThreeVector EndUnitMomentum ;
432 G4double lengthAlongCurve ; 432 G4double lengthAlongCurve ;
433 G4double restMass = pParticleDef->G 433 G4double restMass = pParticleDef->GetPDGMass() ;
434 434
435 fFieldPropagator->SetChargeMomentumMass( 435 fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units
436 momentumMagnitude, // in Mev/c 436 momentumMagnitude, // in Mev/c
437 restMass ) ; 437 restMass ) ;
438 438
439 G4ThreeVector spin = track.GetPola 439 G4ThreeVector spin = track.GetPolarization() ;
440 G4FieldTrack aFieldTrack = G4FieldTrack( 440 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
441 track.GetMomentumDirection(), 441 track.GetMomentumDirection(),
442 0.0, 442 0.0,
443 track.GetKineticEnergy(), 443 track.GetKineticEnergy(),
444 restMass, 444 restMass,
445 track.GetVelocity(), 445 track.GetVelocity(),
446 track.GetGlobalTime(), // Lab. 446 track.GetGlobalTime(), // Lab.
447 track.GetProperTime(), // Part. 447 track.GetProperTime(), // Part.
448 &spin ) ; 448 &spin ) ;
449 if( currentMinimumStep > 0 ) 449 if( currentMinimumStep > 0 )
450 { 450 {
451 // Do the Transport in the field (non rec 451 // Do the Transport in the field (non recti-linear)
452 // 452 //
453 lengthAlongCurve = fFieldPropagator->Comp 453 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
454 currentMinimumStep, 454 currentMinimumStep,
455 currentSafety, 455 currentSafety,
456 track.GetVolume() ) ; 456 track.GetVolume() ) ;
457 State(fGeometryLimitedStep)= lengthAlongC 457 State(fGeometryLimitedStep)= lengthAlongCurve < currentMinimumStep;
458 if( State(fGeometryLimitedStep) ) 458 if( State(fGeometryLimitedStep) )
459 { 459 {
460 geometryStepLength = lengthAlongCurve ; 460 geometryStepLength = lengthAlongCurve ;
461 } 461 }
462 else 462 else
463 { 463 {
464 geometryStepLength = currentMinimumStep 464 geometryStepLength = currentMinimumStep ;
465 } 465 }
466 466
467 // Remember last safety origin & value. 467 // Remember last safety origin & value.
468 // 468 //
469 State(fPreviousSftOrigin) = startPosition 469 State(fPreviousSftOrigin) = startPosition ;
470 State(fPreviousSafety) = currentSafety 470 State(fPreviousSafety) = currentSafety ;
471 fpSafetyHelper->SetCurrentSafety( newSafe 471 fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
472 } 472 }
473 else 473 else
474 { 474 {
475 geometryStepLength = lengthAlongCurve= 475 geometryStepLength = lengthAlongCurve= 0.0 ;
476 State(fGeometryLimitedStep) = false ; 476 State(fGeometryLimitedStep) = false ;
477 } 477 }
478 478
479 // Get the End-Position and End-Momentum 479 // Get the End-Position and End-Momentum (Dir-ection)
480 // 480 //
481 State(fTransportEndPosition) = aFieldTrac 481 State(fTransportEndPosition) = aFieldTrack.GetPosition() ;
482 482
483 // Momentum: Magnitude and direction can 483 // Momentum: Magnitude and direction can be changed too now ...
484 // 484 //
485 State(fMomentumChanged) = true ; 485 State(fMomentumChanged) = true ;
486 State(fTransportEndMomentumDir) = aFieldT 486 State(fTransportEndMomentumDir) = aFieldTrack.GetMomentumDir() ;
487 487
488 State(fTransportEndKineticEnergy) = aFie 488 State(fTransportEndKineticEnergy) = aFieldTrack.GetKineticEnergy() ;
489 489
490 if( fFieldPropagator->GetCurrentFieldMana 490 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
491 { 491 {
492 // If the field can change energy, then t 492 // If the field can change energy, then the time must be integrated
493 // - so this should have been updated 493 // - so this should have been updated
494 // 494 //
495 State(fCandidateEndGlobalTime) = aField 495 State(fCandidateEndGlobalTime) = aFieldTrack.GetLabTimeOfFlight();
496 State(fEndGlobalTimeComputed) = true; 496 State(fEndGlobalTimeComputed) = true;
497 497
498 State(theInteractionTimeLeft) = State(fCa 498 State(theInteractionTimeLeft) = State(fCandidateEndGlobalTime) -
499 track 499 track.GetGlobalTime() ;
500 500
501 // was ( State(fCandidateEndGlobalTime) ! 501 // was ( State(fCandidateEndGlobalTime) != track.GetGlobalTime() );
502 // a cleaner way is to have FieldTrack kn 502 // a cleaner way is to have FieldTrack knowing whether time is updated.
503 } 503 }
504 else 504 else
505 { 505 {
506 // The energy should be unchanged by fiel 506 // The energy should be unchanged by field transport,
507 // - so the time changed will be calcu 507 // - so the time changed will be calculated elsewhere
508 // 508 //
509 State(fEndGlobalTimeComputed) = false; 509 State(fEndGlobalTimeComputed) = false;
510 510
511 // Check that the integration preserved t 511 // Check that the integration preserved the energy
512 // - and if not correct this! 512 // - and if not correct this!
513 G4double startEnergy= track.GetKineticEn 513 G4double startEnergy= track.GetKineticEnergy();
514 G4double endEnergy= State(fTransportEndK 514 G4double endEnergy= State(fTransportEndKineticEnergy);
515 515
516 static G4int no_inexact_steps=0, no_large 516 static G4int no_inexact_steps=0, no_large_ediff;
517 G4double absEdiff = std::fabs(startEnergy 517 G4double absEdiff = std::fabs(startEnergy- endEnergy);
518 if( absEdiff > perMillion * endEnergy ) 518 if( absEdiff > perMillion * endEnergy )
519 { 519 {
520 no_inexact_steps++; 520 no_inexact_steps++;
521 // Possible statistics keeping here ... 521 // Possible statistics keeping here ...
522 } 522 }
523 #ifdef G4VERBOSE 523 #ifdef G4VERBOSE
524 if( fVerboseLevel > 1 ) 524 if( fVerboseLevel > 1 )
525 { 525 {
526 if( std::fabs(startEnergy- endEnergy) > p 526 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
527 { 527 {
528 static G4int no_warnings= 0, warnModulo=1 528 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10;
529 no_large_ediff ++; 529 no_large_ediff ++;
530 if( (no_large_ediff% warnModulo) == 0 ) 530 if( (no_large_ediff% warnModulo) == 0 )
531 { 531 {
532 no_warnings++; 532 no_warnings++;
533 G4cout << "WARNING - G4Transportation::Al 533 G4cout << "WARNING - G4Transportation::AlongStepGetPIL() "
534 << " Energy change in Step is above 1^- 534 << " Energy change in Step is above 1^-3 relative value. " << G4endl
535 << " Relative change in 'tracking' step 535 << " Relative change in 'tracking' step = "
536 << std::setw(15) << (endEnergy-startEnerg 536 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
537 << " Starting E= " << std::setw(12) < 537 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV "
538 << G4endl 538 << G4endl
539 << " Ending E= " << std::setw(12) < 539 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV "
540 << G4endl; 540 << G4endl;
541 G4cout << " Energy has been corrected -- 541 G4cout << " Energy has been corrected -- however, review"
542 << " field propagation parameters for acc 542 << " field propagation parameters for accuracy." << G4endl;
543 if( (fVerboseLevel > 2 ) || (no_warnings< 543 if( (fVerboseLevel > 2 ) || (no_warnings<4) ||
544 (no_large_ediff == warnModulo * modul 544 (no_large_ediff == warnModulo * moduloFactor) )
545 { 545 {
546 G4cout << " These include EpsilonStepMax( 546 G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
547 << " which determine fractional error per 547 << " which determine fractional error per step for integrated quantities. "
548 << G4endl 548 << G4endl
549 << " Note also the influence of the permi 549 << " Note also the influence of the permitted number of integration steps."
550 << G4endl; 550 << G4endl;
551 } 551 }
552 G4cerr << "ERROR - G4Transportation::Alon 552 G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
553 << " Bad 'endpoint'. Energy change 553 << " Bad 'endpoint'. Energy change detected"
554 << " and corrected. " 554 << " and corrected. "
555 << " Has occurred already " 555 << " Has occurred already "
556 << no_large_ediff << " times." << G4endl; 556 << no_large_ediff << " times." << G4endl;
557 if( no_large_ediff == warnModulo * modulo 557 if( no_large_ediff == warnModulo * moduloFactor )
558 { 558 {
559 warnModulo *= moduloFactor; 559 warnModulo *= moduloFactor;
560 } 560 }
561 } 561 }
562 } 562 }
563 } // end of if (fVerboseLevel) 563 } // end of if (fVerboseLevel)
564 #endif 564 #endif
565 // Correct the energy for fields that con 565 // Correct the energy for fields that conserve it
566 // This - hides the integration error 566 // This - hides the integration error
567 // - but gives a better physical an 567 // - but gives a better physical answer
568 State(fTransportEndKineticEnergy)= track. 568 State(fTransportEndKineticEnergy)= track.GetKineticEnergy();
569 } 569 }
570 570
571 State(fTransportEndSpin) = aFieldTrack.Ge 571 State(fTransportEndSpin) = aFieldTrack.GetSpin();
572 State(fParticleIsLooping) = fFieldPropaga 572 State(fParticleIsLooping) = fFieldPropagator->IsParticleLooping() ;
573 State(endpointDistance) = (State(fTrans 573 State(endpointDistance) = (State(fTransportEndPosition) -
574 startPositio 574 startPosition).mag() ;
575 // State(theInteractionTimeLeft) = 575 // State(theInteractionTimeLeft) =
576 track.GetVelocity()/Sta 576 track.GetVelocity()/State(endpointDistance) ;
577 */ 577 */
578 } 578 }
579 579
580 // If we are asked to go a step length of 0, 580 // If we are asked to go a step length of 0, and we are on a boundary
581 // then a boundary will also limit the step 581 // then a boundary will also limit the step -> we must flag this.
582 // 582 //
583 if (currentMinimumStep == 0.0) 583 if (currentMinimumStep == 0.0)
584 { 584 {
585 if (currentSafety == 0.0) 585 if (currentSafety == 0.0)
586 { 586 {
587 State(fGeometryLimitedStep) = true; 587 State(fGeometryLimitedStep) = true;
588 // G4cout << "!!!! Safety is NULL, 588 // G4cout << "!!!! Safety is NULL, on the Boundary !!!!!" << G4endl;
589 // G4cout << " Track position : " < 589 // G4cout << " Track position : " << track.GetPosition() /nanometer
590 // << G4endl; 590 // << G4endl;
591 } 591 }
592 } 592 }
593 593
594 // Update the safety starting from the end-p 594 // Update the safety starting from the end-point,
595 // if it will become negative at the end-poi 595 // if it will become negative at the end-point.
596 // 596 //
597 if (currentSafety < State(fEndPointDistance) 597 if (currentSafety < State(fEndPointDistance))
598 { 598 {
599 // if( particleCharge == 0.0 ) 599 // if( particleCharge == 0.0 )
600 // G4cout << " Avoiding call to Compu 600 // G4cout << " Avoiding call to ComputeSafety : charge = 0.0 " << G4endl;
601 601
602 if (particleCharge != 0.0) 602 if (particleCharge != 0.0)
603 { 603 {
604 604
605 G4double endSafety = fLinearNavigator->C 605 G4double endSafety = fLinearNavigator->ComputeSafety(
606 State(fTransportEndPosition)); 606 State(fTransportEndPosition));
607 currentSafety = endSafety; 607 currentSafety = endSafety;
608 State(fPreviousSftOrigin) = State(fTrans 608 State(fPreviousSftOrigin) = State(fTransportEndPosition);
609 State(fPreviousSafety) = currentSafety; 609 State(fPreviousSafety) = currentSafety;
610 610
611 /* 611 /*
612 G4VTrackStateHandle state = 612 G4VTrackStateHandle state =
613 GetIT(track)->GetTrackingInfo()->GetTra 613 GetIT(track)->GetTrackingInfo()->GetTrackState(fpSafetyHelper);
614 */ 614 */
615 G4TrackStateManager& trackStateMan = Get 615 G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo()
616 ->GetTrackStateManager(); 616 ->GetTrackStateManager();
617 fpSafetyHelper->LoadTrackState(trackStat 617 fpSafetyHelper->LoadTrackState(trackStateMan);
618 // fpSafetyHelper->SetTrackState(state); 618 // fpSafetyHelper->SetTrackState(state);
619 fpSafetyHelper->SetCurrentSafety(current 619 fpSafetyHelper->SetCurrentSafety(currentSafety,
620 State(f 620 State(fTransportEndPosition));
621 fpSafetyHelper->ResetTrackState(); 621 fpSafetyHelper->ResetTrackState();
622 622
623 // Because the Stepping Manager assumes 623 // Because the Stepping Manager assumes it is from the start point,
624 // add the StepLength 624 // add the StepLength
625 // 625 //
626 currentSafety += State(fEndPointDistance 626 currentSafety += State(fEndPointDistance);
627 627
628 #ifdef G4DEBUG_TRANSPORT 628 #ifdef G4DEBUG_TRANSPORT
629 G4cout.precision(12); 629 G4cout.precision(12);
630 G4cout << "***G4Transportation::AlongSte 630 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl;
631 G4cout << " Called Navigator->ComputeSa 631 G4cout << " Called Navigator->ComputeSafety at "
632 << State(fTransportEndPosition) 632 << State(fTransportEndPosition)
633 << " and it returned safety= " << 633 << " and it returned safety= " << endSafety << G4endl;
634 G4cout << " Adding endpoint distance " 634 G4cout << " Adding endpoint distance " << State(fEndPointDistance)
635 << " to obtain pseudo-safety= " 635 << " to obtain pseudo-safety= " << currentSafety << G4endl;
636 #endif 636 #endif
637 } 637 }
638 } 638 }
639 639
640 // fParticleChange.ProposeTrueStepLength(geo 640 // fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
641 641
642 // G4cout << "G4ITTransportation::AlongStepGe 642 // G4cout << "G4ITTransportation::AlongStepGetPhysicalInteractionLength = "
643 // << G4BestUnit(geometryStepLength,"L 643 // << G4BestUnit(geometryStepLength,"Length") << G4endl;
644 644
645 return geometryStepLength; 645 return geometryStepLength;
646 } 646 }
647 647
648 void G4ITTransportation::ComputeStep(const G4T 648 void G4ITTransportation::ComputeStep(const G4Track& track,
649 const G4S 649 const G4Step& /*step*/,
650 const dou 650 const double timeStep,
651 double& o 651 double& oPhysicalStep)
652 { 652 {
653 PrepareState(); << 653 PrepareState()
654 const G4DynamicParticle* pParticle = track.G 654 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
655 G4ThreeVector startMomentumDir = pParticle-> 655 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection();
656 G4ThreeVector startPosition = track.GetPosit 656 G4ThreeVector startPosition = track.GetPosition();
657 657
658 track.CalculateVelocity(); 658 track.CalculateVelocity();
659 G4double initialVelocity = track.GetVelocity 659 G4double initialVelocity = track.GetVelocity();
660 660
661 State(fGeometryLimitedStep) = false; 661 State(fGeometryLimitedStep) = false;
662 662
663 ///////////////////////// 663 /////////////////////////
664 // !!! CASE NO FIELD !!! 664 // !!! CASE NO FIELD !!!
665 ///////////////////////// 665 /////////////////////////
666 State(fCandidateEndGlobalTime) = timeStep + 666 State(fCandidateEndGlobalTime) = timeStep + track.GetGlobalTime();
667 State(fEndGlobalTimeComputed) = true; 667 State(fEndGlobalTimeComputed) = true;
668 668
669 // Choose the calculation of the transportat 669 // Choose the calculation of the transportation: Field or not
670 // 670 //
671 if (!State(fMomentumChanged)) 671 if (!State(fMomentumChanged))
672 { 672 {
673 // G4cout << "Momentum has not chan 673 // G4cout << "Momentum has not changed" << G4endl;
674 fParticleChange.ProposeVelocity(initialVel 674 fParticleChange.ProposeVelocity(initialVelocity);
675 oPhysicalStep = initialVelocity * timeStep 675 oPhysicalStep = initialVelocity * timeStep;
676 676
677 // Calculate final position 677 // Calculate final position
678 // 678 //
679 State(fTransportEndPosition) = startPositi 679 State(fTransportEndPosition) = startPosition
680 + oPhysicalStep * startMomentumDir; 680 + oPhysicalStep * startMomentumDir;
681 } 681 }
682 } 682 }
683 683
684 ////////////////////////////////////////////// 684 //////////////////////////////////////////////////////////////////////////
685 // 685 //
686 // Initialize ParticleChange (by setting al 686 // Initialize ParticleChange (by setting all its members equal
687 // to correspond 687 // to corresponding members in G4Track)
688 #include "G4ParticleTable.hh" 688 #include "G4ParticleTable.hh"
689 G4VParticleChange* G4ITTransportation::AlongSt 689 G4VParticleChange* G4ITTransportation::AlongStepDoIt(const G4Track& track,
690 690 const G4Step& stepData)
691 { 691 {
692 692
693 #if defined (DEBUG_MEM) 693 #if defined (DEBUG_MEM)
694 MemStat mem_first, mem_second, mem_diff; 694 MemStat mem_first, mem_second, mem_diff;
695 #endif 695 #endif
696 696
697 #if defined (DEBUG_MEM) 697 #if defined (DEBUG_MEM)
698 mem_first = MemoryUsage(); 698 mem_first = MemoryUsage();
699 #endif 699 #endif
700 700
701 PrepareState(); << 701 PrepareState()
702 << 702
703 // G4cout << "G4ITTransportation::AlongStepD 703 // G4cout << "G4ITTransportation::AlongStepDoIt" << G4endl;
704 // set pdefOpticalPhoton 704 // set pdefOpticalPhoton
705 // Andrea Dotti: the following statement sho 705 // Andrea Dotti: the following statement should be in a single line:
706 // G4-MT transformation tools get confused i 706 // G4-MT transformation tools get confused if statement spans two lines
707 // If needed contact: adotti@slac.stanford.e 707 // If needed contact: adotti@slac.stanford.edu
708 static G4ThreadLocal G4ParticleDefinition* p << 708 static G4ThreadLocal G4ParticleDefinition* pdefOpticalPhoton = 0;
709 if (pdefOpticalPhoton == nullptr) pdefOptica << 709 if (!pdefOpticalPhoton) pdefOpticalPhoton =
710 G4ParticleTable::GetParticleTable()->Fin 710 G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
711 711
712 static G4ThreadLocal G4int noCalls = 0; 712 static G4ThreadLocal G4int noCalls = 0;
713 noCalls++; 713 noCalls++;
714 714
715 fParticleChange.Initialize(track); 715 fParticleChange.Initialize(track);
716 716
717 // Code for specific process 717 // Code for specific process
718 // 718 //
719 fParticleChange.ProposePosition(State(fTrans 719 fParticleChange.ProposePosition(State(fTransportEndPosition));
720 fParticleChange.ProposeMomentumDirection(Sta 720 fParticleChange.ProposeMomentumDirection(State(fTransportEndMomentumDir));
721 fParticleChange.ProposeEnergy(State(fTranspo 721 fParticleChange.ProposeEnergy(State(fTransportEndKineticEnergy));
722 fParticleChange.SetMomentumChanged(State(fMo 722 fParticleChange.SetMomentumChanged(State(fMomentumChanged));
723 723
724 fParticleChange.ProposePolarization(State(fT 724 fParticleChange.ProposePolarization(State(fTransportEndSpin));
725 725
726 G4double deltaTime = 0.0; 726 G4double deltaTime = 0.0;
727 727
728 // Calculate Lab Time of Flight (ONLY if fi 728 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
729 // G4double endTime = State(fCandidateEndG 729 // G4double endTime = State(fCandidateEndGlobalTime);
730 // G4double delta_time = endTime - startTime 730 // G4double delta_time = endTime - startTime;
731 731
732 G4double startTime = track.GetGlobalTime(); 732 G4double startTime = track.GetGlobalTime();
733 ///_________________________________________ 733 ///___________________________________________________________________________
734 /// !!!!!!! 734 /// !!!!!!!
735 /// A REVOIR !!!! 735 /// A REVOIR !!!!
736 if (State(fEndGlobalTimeComputed) == 0) << 736 if (State(fEndGlobalTimeComputed) == false)
737 { 737 {
738 // The time was not integrated .. make the 738 // The time was not integrated .. make the best estimate possible
739 // 739 //
740 G4double initialVelocity = stepData.GetPre 740 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
741 G4double stepLength = track.GetStepLength( 741 G4double stepLength = track.GetStepLength();
742 742
743 deltaTime = 0.0; // in case initialVelocit 743 deltaTime = 0.0; // in case initialVelocity = 0
744 if (track.GetParticleDefinition() == pdefO 744 if (track.GetParticleDefinition() == pdefOpticalPhoton)
745 { 745 {
746 // For only Optical Photon, final veloci 746 // For only Optical Photon, final velocity is used
747 double finalVelocity = track.CalculateVe 747 double finalVelocity = track.CalculateVelocityForOpticalPhoton();
748 fParticleChange.ProposeVelocity(finalVel 748 fParticleChange.ProposeVelocity(finalVelocity);
749 deltaTime = stepLength / finalVelocity; 749 deltaTime = stepLength / finalVelocity;
750 } 750 }
751 else if (initialVelocity > 0.0) 751 else if (initialVelocity > 0.0)
752 { 752 {
753 deltaTime = stepLength / initialVelocity 753 deltaTime = stepLength / initialVelocity;
754 } 754 }
755 755
756 State(fCandidateEndGlobalTime) = startTime 756 State(fCandidateEndGlobalTime) = startTime + deltaTime;
757 } 757 }
758 else 758 else
759 { 759 {
760 deltaTime = State(fCandidateEndGlobalTime) 760 deltaTime = State(fCandidateEndGlobalTime) - startTime;
761 } 761 }
762 762
763 fParticleChange.ProposeGlobalTime(State(fCan 763 fParticleChange.ProposeGlobalTime(State(fCandidateEndGlobalTime));
764 fParticleChange.ProposeLocalTime(track.GetLo 764 fParticleChange.ProposeLocalTime(track.GetLocalTime() + deltaTime);
765 /* 765 /*
766 // Now Correct by Lorentz factor to get del 766 // Now Correct by Lorentz factor to get delta "proper" Time
767 767
768 G4double restMass = track.GetDynamic 768 G4double restMass = track.GetDynamicParticle()->GetMass() ;
769 G4double deltaProperTime = deltaTime*( rest 769 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
770 770
771 fParticleChange.ProposeProperTime(track.Get 771 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
772 */ 772 */
773 773
774 fParticleChange.ProposeTrueStepLength(track. 774 fParticleChange.ProposeTrueStepLength(track.GetStepLength());
775 775
776 ///_________________________________________ 776 ///___________________________________________________________________________
777 /// 777 ///
778 778
779 // If the particle is caught looping or is s 779 // If the particle is caught looping or is stuck (in very difficult
780 // boundaries) in a magnetic field (doing ma 780 // boundaries) in a magnetic field (doing many steps)
781 // THEN this kills it ... 781 // THEN this kills it ...
782 // 782 //
783 if (State(fParticleIsLooping)) 783 if (State(fParticleIsLooping))
784 { 784 {
785 G4double endEnergy = State(fTransportEndKi 785 G4double endEnergy = State(fTransportEndKineticEnergy);
786 786
787 if ((endEnergy < fThreshold_Important_Ener 787 if ((endEnergy < fThreshold_Important_Energy) || (State(fNoLooperTrials)
788 >= fThresholdTrials)) 788 >= fThresholdTrials))
789 { 789 {
790 // Kill the looping particle 790 // Kill the looping particle
791 // 791 //
792 // G4cout << "G4ITTransportation will ki 792 // G4cout << "G4ITTransportation will killed the molecule"<< G4endl;
793 fParticleChange.ProposeTrackStatus(fStop 793 fParticleChange.ProposeTrackStatus(fStopAndKill);
794 794
795 // 'Bare' statistics 795 // 'Bare' statistics
796 fSumEnergyKilled += endEnergy; 796 fSumEnergyKilled += endEnergy;
797 if (endEnergy > fMaxEnergyKilled) 797 if (endEnergy > fMaxEnergyKilled)
798 { 798 {
799 fMaxEnergyKilled = endEnergy; 799 fMaxEnergyKilled = endEnergy;
800 } 800 }
801 801
802 #ifdef G4VERBOSE 802 #ifdef G4VERBOSE
803 if ((fVerboseLevel > 1) || (endEnergy > 803 if ((fVerboseLevel > 1) || (endEnergy > fThreshold_Warning_Energy))
804 { 804 {
805 G4cout 805 G4cout
806 << " G4ITTransportation is killing 806 << " G4ITTransportation is killing track that is looping or stuck "
807 << G4endl<< " This track has " < 807 << G4endl<< " This track has " << track.GetKineticEnergy() / MeV
808 << " MeV energy." << G4endl; 808 << " MeV energy." << G4endl;
809 G4cout << " Number of trials = " << 809 G4cout << " Number of trials = " << State(fNoLooperTrials)
810 << " No of calls to AlongStepDoIt = 810 << " No of calls to AlongStepDoIt = " << noCalls
811 << G4endl; 811 << G4endl;
812 } 812 }
813 #endif 813 #endif
814 State(fNoLooperTrials) = 0; 814 State(fNoLooperTrials) = 0;
815 } 815 }
816 else 816 else
817 { 817 {
818 State(fNoLooperTrials)++; 818 State(fNoLooperTrials)++;
819 #ifdef G4VERBOSE 819 #ifdef G4VERBOSE
820 if ((fVerboseLevel > 2)) 820 if ((fVerboseLevel > 2))
821 { 821 {
822 G4cout << " G4ITTransportation::Alon 822 G4cout << " G4ITTransportation::AlongStepDoIt(): Particle looping - "
823 << " Number of trials = " << 823 << " Number of trials = " << State(fNoLooperTrials)
824 << " No of calls to = " << n 824 << " No of calls to = " << noCalls << G4endl;
825 } 825 }
826 #endif 826 #endif
827 } 827 }
828 } 828 }
829 else 829 else
830 { 830 {
831 State(fNoLooperTrials)=0; 831 State(fNoLooperTrials)=0;
832 } 832 }
833 833
834 // Another (sometimes better way) is to use 834 // Another (sometimes better way) is to use a user-limit maximum Step size
835 // to alleviate this problem .. 835 // to alleviate this problem ..
836 836
837 // Introduce smooth curved trajectories to p 837 // Introduce smooth curved trajectories to particle-change
838 // 838 //
839 fParticleChange.SetPointerToVectorOfAuxiliar 839 fParticleChange.SetPointerToVectorOfAuxiliaryPoints(
840 fFieldPropagator->GimmeTrajectoryVectorA 840 fFieldPropagator->GimmeTrajectoryVectorAndForgetIt());
841 841
842 #if defined (DEBUG_MEM) 842 #if defined (DEBUG_MEM)
843 mem_second = MemoryUsage(); 843 mem_second = MemoryUsage();
844 mem_diff = mem_second-mem_first; 844 mem_diff = mem_second-mem_first;
845 G4cout << "\t || MEM || End of G4ITTransport 845 G4cout << "\t || MEM || End of G4ITTransportation::AlongStepDoIt, diff is: "
846 << mem_diff << G4endl; 846 << mem_diff << G4endl;
847 #endif 847 #endif
848 848
849 return &fParticleChange; 849 return &fParticleChange;
850 } 850 }
851 851
852 ////////////////////////////////////////////// 852 //////////////////////////////////////////////////////////////////////////
853 // 853 //
854 // This ensures that the PostStep action is a 854 // This ensures that the PostStep action is always called,
855 // so that it can do the relocation if it is 855 // so that it can do the relocation if it is needed.
856 // 856 //
857 857
858 G4double 858 G4double
859 G4ITTransportation:: 859 G4ITTransportation::
860 PostStepGetPhysicalInteractionLength(const G4T 860 PostStepGetPhysicalInteractionLength(const G4Track&, // track
861 G4double, 861 G4double, // previousStepSize
862 G4ForceCo 862 G4ForceCondition* pForceCond)
863 { 863 {
864 *pForceCond = Forced; 864 *pForceCond = Forced;
865 return DBL_MAX; // was kInfinity ; but conve 865 return DBL_MAX; // was kInfinity ; but convention now is DBL_MAX
866 } 866 }
867 867
868 ////////////////////////////////////////////// 868 /////////////////////////////////////////////////////////////////////////////
869 // 869 //
870 870
871 G4VParticleChange* G4ITTransportation::PostSte 871 G4VParticleChange* G4ITTransportation::PostStepDoIt(const G4Track& track,
872 872 const G4Step&)
873 { 873 {
874 // G4cout << "G4ITTransportation::PostSte 874 // G4cout << "G4ITTransportation::PostStepDoIt" << G4endl;
875 << 875
876 PrepareState(); << 876 PrepareState()
877 G4TouchableHandle retCurrentTouchable; // Th 877 G4TouchableHandle retCurrentTouchable; // The one to return
878 G4bool isLastStep = false; 878 G4bool isLastStep = false;
879 879
880 // Initialize ParticleChange (by setting al 880 // Initialize ParticleChange (by setting all its members equal
881 // to correspond 881 // to corresponding members in G4Track)
882 fParticleChange.Initialize(track); // To ini 882 fParticleChange.Initialize(track); // To initialise TouchableChange
883 883
884 fParticleChange.ProposeTrackStatus(track.Get 884 fParticleChange.ProposeTrackStatus(track.GetTrackStatus());
885 885
886 // If the Step was determined by the volume 886 // If the Step was determined by the volume boundary,
887 // logically relocate the particle 887 // logically relocate the particle
888 888
889 if (State(fGeometryLimitedStep)) 889 if (State(fGeometryLimitedStep))
890 { 890 {
891 891
892 if(fVerboseLevel != 0) << 892 if(fVerboseLevel)
893 { 893 {
894 G4cout << "Step is limited by geometry " 894 G4cout << "Step is limited by geometry "
895 << "track ID : " << track.GetTrac 895 << "track ID : " << track.GetTrackID() << G4endl;
896 } 896 }
897 897
898 // fCurrentTouchable will now become the p 898 // fCurrentTouchable will now become the previous touchable,
899 // and what was the previous will be freed 899 // and what was the previous will be freed.
900 // (Needed because the preStepPoint can po 900 // (Needed because the preStepPoint can point to the previous touchable)
901 901
902 if ( State(fCurrentTouchableHandle)->GetVo << 902 if ( State(fCurrentTouchableHandle)->GetVolume() == 0)
903 { 903 {
904 G4ExceptionDescription exceptionDescript 904 G4ExceptionDescription exceptionDescription;
905 exceptionDescription << "No current touc 905 exceptionDescription << "No current touchable found ";
906 G4Exception(" G4ITTransportation::PostSt 906 G4Exception(" G4ITTransportation::PostStepDoIt", "G4ITTransportation001",
907 FatalErrorInArgument, except 907 FatalErrorInArgument, exceptionDescription);
908 } 908 }
909 909
910 fLinearNavigator->SetGeometricallyLimitedS 910 fLinearNavigator->SetGeometricallyLimitedStep();
911 fLinearNavigator->LocateGlobalPointAndUpda 911 fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
912 track.GetPosition(), track.GetMomentum 912 track.GetPosition(), track.GetMomentumDirection(),
913 State(fCurrentTouchableHandle), true); 913 State(fCurrentTouchableHandle), true);
914 // Check whether the particle is out of th 914 // Check whether the particle is out of the world volume
915 // If so it has exited and must be killed. 915 // If so it has exited and must be killed.
916 // 916 //
917 if ( State(fCurrentTouchableHandle)->GetVo << 917 if ( State(fCurrentTouchableHandle)->GetVolume() == 0)
918 { 918 {
919 // abort(); 919 // abort();
920 #ifdef G4VERBOSE 920 #ifdef G4VERBOSE
921 if (fVerboseLevel > 0) 921 if (fVerboseLevel > 0)
922 { 922 {
923 G4cout << "Track position : " << track 923 G4cout << "Track position : " << track.GetPosition() / nanometer
924 << " [nm]" << " Track ID : " << 924 << " [nm]" << " Track ID : " << track.GetTrackID() << G4endl;
925 G4cout << "G4ITTransportation will kil 925 G4cout << "G4ITTransportation will killed the track because "
926 "State(fCurrentTouchableHandle)->G 926 "State(fCurrentTouchableHandle)->GetVolume() == 0"<< G4endl;
927 } 927 }
928 #endif 928 #endif
929 fParticleChange.ProposeTrackStatus( fSto 929 fParticleChange.ProposeTrackStatus( fStopAndKill );
930 } 930 }
931 931
932 retCurrentTouchable = State(fCurrentToucha 932 retCurrentTouchable = State(fCurrentTouchableHandle);
933 933
934 // G4cout << "Current volume : " << tra 934 // G4cout << "Current volume : " << track.GetVolume()->GetName()
935 // << " Next volume : " 935 // << " Next volume : "
936 // << (State(fCurrentTouchableHa 936 // << (State(fCurrentTouchableHandle)->GetVolume() ?
937 // State(fCurrentTouchableHandle)->GetVolum 937 // State(fCurrentTouchableHandle)->GetVolume()->GetName():"OutWorld")
938 // << " Position : " << track.Ge 938 // << " Position : " << track.GetPosition() / nanometer
939 // << " track ID : " << track.Ge 939 // << " track ID : " << track.GetTrackID()
940 // << G4endl; 940 // << G4endl;
941 941
942 fParticleChange.SetTouchableHandle(State(f 942 fParticleChange.SetTouchableHandle(State(fCurrentTouchableHandle));
943 943
944 // Update the Step flag which identifies t 944 // Update the Step flag which identifies the Last Step in a volume
945 isLastStep = fLinearNavigator->ExitedMothe 945 isLastStep = fLinearNavigator->ExitedMotherVolume()
946 || fLinearNavigator->EnteredDau << 946 | fLinearNavigator->EnteredDaughterVolume();
947 947
948 #ifdef G4DEBUG_TRANSPORT 948 #ifdef G4DEBUG_TRANSPORT
949 // Checking first implementation of flagg 949 // Checking first implementation of flagging Last Step in Volume
950 G4bool exiting = fLinearNavigator->ExitedM 950 G4bool exiting = fLinearNavigator->ExitedMotherVolume();
951 G4bool entering = fLinearNavigator->Entere 951 G4bool entering = fLinearNavigator->EnteredDaughterVolume();
952 952
953 if( ! (exiting || entering) ) 953 if( ! (exiting || entering) )
954 { 954 {
955 G4cout << " Transport> : Proposed isLas 955 G4cout << " Transport> : Proposed isLastStep= " << isLastStep
956 << " Exiting " << fLinearNavigator->Exit 956 << " Exiting " << fLinearNavigator->ExitedMotherVolume()
957 << " Entering " << fLinearNavigator->Ent 957 << " Entering " << fLinearNavigator->EnteredDaughterVolume()
958 << " Track position : " << track.GetPosi 958 << " Track position : " << track.GetPosition() /nanometer << " [nm]"
959 << G4endl; 959 << G4endl;
960 G4cout << " Track position : " << track. 960 G4cout << " Track position : " << track.GetPosition() /nanometer
961 << G4endl; 961 << G4endl;
962 } 962 }
963 #endif 963 #endif
964 } 964 }
965 else // fGeometryLimitedStep is false 965 else // fGeometryLimitedStep is false
966 { 966 {
967 // This serves only to move the Navigator' 967 // This serves only to move the Navigator's location
968 // 968 //
969 // abort(); 969 // abort();
970 fLinearNavigator->LocateGlobalPointWithinV 970 fLinearNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
971 971
972 // The value of the track's current Toucha 972 // The value of the track's current Touchable is retained.
973 // (and it must be correct because we must 973 // (and it must be correct because we must use it below to
974 // overwrite the (unset) one in particle c 974 // overwrite the (unset) one in particle change)
975 // It must be fCurrentTouchable too ?? 975 // It must be fCurrentTouchable too ??
976 // 976 //
977 fParticleChange.SetTouchableHandle(track.G 977 fParticleChange.SetTouchableHandle(track.GetTouchableHandle());
978 retCurrentTouchable = track.GetTouchableHa 978 retCurrentTouchable = track.GetTouchableHandle();
979 979
980 isLastStep = false; 980 isLastStep = false;
981 #ifdef G4DEBUG_TRANSPORT 981 #ifdef G4DEBUG_TRANSPORT
982 // Checking first implementation of flagg 982 // Checking first implementation of flagging Last Step in Volume
983 G4cout << " Transport> Proposed isLastStep 983 G4cout << " Transport> Proposed isLastStep= " << isLastStep
984 << " Geometry did not limit step. Position 984 << " Geometry did not limit step. Position : "
985 << track.GetPosition()/ nanometer << G4end 985 << track.GetPosition()/ nanometer << G4endl;
986 #endif 986 #endif
987 } // endif ( fGeometryLimitedStep ) 987 } // endif ( fGeometryLimitedStep )
988 988
989 fParticleChange.ProposeLastStepInVolume(isLa 989 fParticleChange.ProposeLastStepInVolume(isLastStep);
990 990
991 const G4VPhysicalVolume* pNewVol = retCurren 991 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume();
992 const G4Material* pNewMaterial = nullptr; << 992 const G4Material* pNewMaterial = 0;
993 G4VSensitiveDetector* pNewSensitiveDetector << 993 const G4VSensitiveDetector* pNewSensitiveDetector = 0;
994 994
995 if (pNewVol != nullptr) << 995 if (pNewVol != 0)
996 { 996 {
997 pNewMaterial = pNewVol->GetLogicalVolume() 997 pNewMaterial = pNewVol->GetLogicalVolume()->GetMaterial();
998 pNewSensitiveDetector = pNewVol->GetLogica 998 pNewSensitiveDetector = pNewVol->GetLogicalVolume()->GetSensitiveDetector();
999 } 999 }
1000 1000
1001 // ( <const_cast> pNewMaterial ) ; 1001 // ( <const_cast> pNewMaterial ) ;
>> 1002 // ( <const_cast> pNewSensitiveDetector) ;
1002 1003
1003 fParticleChange.SetMaterialInTouchable((G4M 1004 fParticleChange.SetMaterialInTouchable((G4Material *) pNewMaterial);
1004 fParticleChange.SetSensitiveDetectorInTouch << 1005 fParticleChange.SetSensitiveDetectorInTouchable(
>> 1006 (G4VSensitiveDetector *) pNewSensitiveDetector);
1005 1007
1006 const G4MaterialCutsCouple* pNewMaterialCut << 1008 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
1007 if (pNewVol != nullptr) << 1009 if (pNewVol != 0)
1008 { 1010 {
1009 pNewMaterialCutsCouple = 1011 pNewMaterialCutsCouple =
1010 pNewVol->GetLogicalVolume()->GetMater 1012 pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
1011 } 1013 }
1012 1014
1013 if (pNewVol != nullptr && pNewMaterialCutsC << 1015 if (pNewVol != 0 && pNewMaterialCutsCouple != 0
1014 && pNewMaterialCutsCouple->GetMaterial( 1016 && pNewMaterialCutsCouple->GetMaterial() != pNewMaterial)
1015 { 1017 {
1016 // for parametrized volume 1018 // for parametrized volume
1017 // 1019 //
1018 pNewMaterialCutsCouple = G4ProductionCuts 1020 pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()
1019 ->GetMaterialCutsCouple(pNewMaterial, 1021 ->GetMaterialCutsCouple(pNewMaterial,
1020 pNewMaterialC 1022 pNewMaterialCutsCouple->GetProductionCuts());
1021 } 1023 }
1022 fParticleChange.SetMaterialCutsCoupleInTouc 1024 fParticleChange.SetMaterialCutsCoupleInTouchable(pNewMaterialCutsCouple);
1023 1025
1024 // temporarily until Get/Set Material of Pa 1026 // temporarily until Get/Set Material of ParticleChange,
1025 // and StepPoint can be made const. 1027 // and StepPoint can be made const.
1026 // Set the touchable in ParticleChange 1028 // Set the touchable in ParticleChange
1027 // this must always be done because the par 1029 // this must always be done because the particle change always
1028 // uses this value to overwrite the current 1030 // uses this value to overwrite the current touchable pointer.
1029 // 1031 //
1030 fParticleChange.SetTouchableHandle(retCurre 1032 fParticleChange.SetTouchableHandle(retCurrentTouchable);
1031 1033
1032 return &fParticleChange; 1034 return &fParticleChange;
1033 } 1035 }
1034 1036
1035 // New method takes over the responsibility t 1037 // New method takes over the responsibility to reset the state of
1036 // G4Transportation object at the start of a 1038 // G4Transportation object at the start of a new track or the resumption of
1037 // a suspended track. 1039 // a suspended track.
1038 1040
1039 void G4ITTransportation::StartTracking(G4Trac 1041 void G4ITTransportation::StartTracking(G4Track* track)
1040 { 1042 {
1041 G4VProcess::StartTracking(track); 1043 G4VProcess::StartTracking(track);
1042 if (fInstantiateProcessState) 1044 if (fInstantiateProcessState)
1043 { 1045 {
1044 // G4VITProcess::fpState = new G4ITTra 1046 // G4VITProcess::fpState = new G4ITTransportationState();
1045 G4VITProcess::fpState = std::make_shared< << 1047 G4VITProcess::fpState.reset(new G4ITTransportationState());
1046 // Will set in the same time fTransportat 1048 // Will set in the same time fTransportationState
1047 } 1049 }
1048 1050
1049 fpSafetyHelper->NewTrackState(); 1051 fpSafetyHelper->NewTrackState();
1050 fpSafetyHelper->SaveTrackState( 1052 fpSafetyHelper->SaveTrackState(
1051 GetIT(track)->GetTrackingInfo()->GetTra 1053 GetIT(track)->GetTrackingInfo()->GetTrackStateManager());
1052 1054
1053 // The actions here are those that were tak 1055 // The actions here are those that were taken in AlongStepGPIL
1054 // when track.GetCurrentStepNumber()==1 1056 // when track.GetCurrentStepNumber()==1
1055 1057
1056 // reset safety value and center 1058 // reset safety value and center
1057 // 1059 //
1058 // State(fPreviousSafety) = 0.0 ; 1060 // State(fPreviousSafety) = 0.0 ;
1059 // State(fPreviousSftOrigin) = G4ThreeVe 1061 // State(fPreviousSftOrigin) = G4ThreeVector(0.,0.,0.) ;
1060 1062
1061 // reset looping counter -- for motion in f 1063 // reset looping counter -- for motion in field
1062 // State(fNoLooperTrials)= 0; 1064 // State(fNoLooperTrials)= 0;
1063 // Must clear this state .. else it depends 1065 // Must clear this state .. else it depends on last track's value
1064 // --> a better solution would set this fr 1066 // --> a better solution would set this from state of suspended track TODO ?
1065 // Was if( aTrack->GetCurrentStepNumber()== 1067 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
1066 1068
1067 // ChordFinder reset internal state 1069 // ChordFinder reset internal state
1068 // 1070 //
1069 if (DoesGlobalFieldExist()) 1071 if (DoesGlobalFieldExist())
1070 { 1072 {
1071 fFieldPropagator->ClearPropagatorState(); 1073 fFieldPropagator->ClearPropagatorState();
1072 // Resets all state of field propagator c 1074 // Resets all state of field propagator class (ONLY)
1073 // including safety values (in case of ov 1075 // including safety values (in case of overlaps and to wipe for first track).
1074 1076
1075 // G4ChordFinder* chordF= fFieldPropagato 1077 // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
1076 // if( chordF ) chordF->ResetStepEstimate 1078 // if( chordF ) chordF->ResetStepEstimate();
1077 } 1079 }
1078 1080
1079 // Make sure to clear the chord finders of 1081 // Make sure to clear the chord finders of all fields (ie managers)
1080 static G4ThreadLocal G4FieldManagerStore* f << 1082 static G4ThreadLocal G4FieldManagerStore* fieldMgrStore = 0;
1081 if (fieldMgrStore == nullptr) fieldMgrStore << 1083 if (!fieldMgrStore) fieldMgrStore = G4FieldManagerStore::GetInstance();
1082 fieldMgrStore->ClearAllChordFindersState(); 1084 fieldMgrStore->ClearAllChordFindersState();
1083 1085
1084 // Update the current touchable handle (fr 1086 // Update the current touchable handle (from the track's)
1085 // 1087 //
1086 PrepareState(); << 1088 PrepareState()
1087 State(fCurrentTouchableHandle) = track->Get 1089 State(fCurrentTouchableHandle) = track->GetTouchableHandle();
1088 1090
1089 G4VITProcess::StartTracking(track); 1091 G4VITProcess::StartTracking(track);
1090 } 1092 }
1091 1093
1092 #undef State 1094 #undef State
1093 #undef PrepareState 1095 #undef PrepareState
1094 1096