| 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 // G4BiasingProcessInterface <<
27 // ------------------------------------------- <<
28 <<
29 #include "G4BiasingProcessInterface.hh" 26 #include "G4BiasingProcessInterface.hh"
30 #include "G4VBiasingOperator.hh" 27 #include "G4VBiasingOperator.hh"
31 #include "G4VBiasingOperation.hh" 28 #include "G4VBiasingOperation.hh"
32 #include "G4ParticleChangeForOccurenceBiasing. 29 #include "G4ParticleChangeForOccurenceBiasing.hh"
33 #include "G4ParticleChangeForNothing.hh" 30 #include "G4ParticleChangeForNothing.hh"
34 #include "G4VBiasingInteractionLaw.hh" 31 #include "G4VBiasingInteractionLaw.hh"
35 #include "G4InteractionLawPhysical.hh" 32 #include "G4InteractionLawPhysical.hh"
36 #include "G4ProcessManager.hh" 33 #include "G4ProcessManager.hh"
37 #include "G4BiasingAppliedCase.hh" 34 #include "G4BiasingAppliedCase.hh"
38 #include "G4ParallelGeometriesLimiterProcess.h 35 #include "G4ParallelGeometriesLimiterProcess.hh"
39 36
40 G4Cache<G4bool> G4BiasingProcessInterface::fRe << 37 G4Cache<G4bool> G4BiasingProcessInterface::fResetInteractionLaws;// = true;
41 G4Cache<G4bool> G4BiasingProcessInterface::fCo << 38 G4Cache<G4bool> G4BiasingProcessInterface::fCommonStart;// = true;
42 G4Cache<G4bool> G4BiasingProcessInterface::fCo << 39 G4Cache<G4bool> G4BiasingProcessInterface::fCommonEnd;// = true;
43 G4Cache<G4bool> G4BiasingProcessInterface::fDo << 40 G4Cache<G4bool> G4BiasingProcessInterface::fDoCommonConfigure;
44 << 41
45 G4BiasingProcessInterface::G4BiasingProcessInt << 42 G4BiasingProcessInterface::G4BiasingProcessInterface(G4String name)
46 : G4VProcess( name ), << 43 : G4VProcess ( name ),
47 fResetWrappedProcessInteractionLength( tru << 44 fCurrentTrack ( nullptr ),
>> 45 fPreviousStepSize (-1.0), fCurrentMinimumStep( -1.0 ), fProposedSafety ( -1.0),
>> 46 fOccurenceBiasingOperation( nullptr ), fFinalStateBiasingOperation( nullptr ), fNonPhysicsBiasingOperation( nullptr ),
>> 47 fPreviousOccurenceBiasingOperation( nullptr ), fPreviousFinalStateBiasingOperation( nullptr ), fPreviousNonPhysicsBiasingOperation( nullptr ),
>> 48 fResetWrappedProcessInteractionLength( true ),
>> 49 fWrappedProcess ( nullptr ),
>> 50 fIsPhysicsBasedBiasing ( false ),
>> 51 fWrappedProcessIsAtRest ( false ),
>> 52 fWrappedProcessIsAlong ( false ),
>> 53 fWrappedProcessIsPost ( false ),
>> 54 fWrappedProcessPostStepGPIL ( -1.0 ),
>> 55 fBiasingPostStepGPIL ( -1.0 ),
>> 56 fWrappedProcessInteractionLength ( -1.0 ),
>> 57 fWrappedProcessForceCondition ( NotForced ),
>> 58 fBiasingForceCondition ( NotForced ),
>> 59 fWrappedProcessAlongStepGPIL ( -1.0 ),
>> 60 fBiasingAlongStepGPIL ( -1.0 ),
>> 61 fWrappedProcessGPILSelection ( NotCandidateForSelection ),
>> 62 fBiasingGPILSelection ( NotCandidateForSelection ),
>> 63 fBiasingInteractionLaw ( nullptr ),
>> 64 fPreviousBiasingInteractionLaw ( nullptr ),
>> 65 fPhysicalInteractionLaw ( nullptr ),
>> 66 fOccurenceBiasingParticleChange ( nullptr ),
>> 67 fDummyParticleChange ( nullptr ),
>> 68 fIamFirstGPIL ( false ),
>> 69 fProcessManager ( nullptr ),
>> 70 fSharedData ( nullptr )
48 { 71 {
49 for (G4int i = 0 ; i < 8 ; i++) fFirstLastF 72 for (G4int i = 0 ; i < 8 ; i++) fFirstLastFlags[i] = false;
50 fResetInteractionLaws.Put( true ); 73 fResetInteractionLaws.Put( true );
51 fCommonStart .Put( true ); << 74 fCommonStart.Put(true);
52 fCommonEnd .Put( true ); << 75 fCommonEnd.Put(true);
53 fDoCommonConfigure .Put( true ); << 76 fDoCommonConfigure.Put(true);
54 } 77 }
55 78
56 79
57 G4BiasingProcessInterface:: << 80 G4BiasingProcessInterface::G4BiasingProcessInterface(G4VProcess* wrappedProcess,
58 G4BiasingProcessInterface(G4VProcess* wrappedP << 81 G4bool wrappedIsAtRest, G4bool wrappedIsAlongStep, G4bool wrappedIsPostStep,
59 G4bool wrappedIsAtRe << 82 G4String useThisName)
60 G4bool wrappedIsAlon << 83 : G4VProcess( useThisName != "" ? useThisName : "biasWrapper("+wrappedProcess->GetProcessName()+")",
61 G4bool wrappedIsPost << 84 wrappedProcess->GetProcessType()),
62 G4String useThisName << 85 fCurrentTrack ( nullptr ),
63 : G4VProcess(useThisName != "" << 86 fPreviousStepSize (-1.0), fCurrentMinimumStep( -1.0 ), fProposedSafety ( -1.0),
64 ? useThisName << 87 fOccurenceBiasingOperation( nullptr ), fFinalStateBiasingOperation( nullptr ), fNonPhysicsBiasingOperation( nullptr ),
65 : "biasWrapper("+wrappedProce << 88 fPreviousOccurenceBiasingOperation( nullptr ), fPreviousFinalStateBiasingOperation( nullptr ), fPreviousNonPhysicsBiasingOperation( nullptr ),
66 wrappedProcess->GetProcessType( << 89 fResetWrappedProcessInteractionLength ( false ),
67 fWrappedProcess( wrappedProcess ), << 90 fWrappedProcess ( wrappedProcess ),
68 fIsPhysicsBasedBiasing( true ), << 91 fIsPhysicsBasedBiasing ( true ),
69 fWrappedProcessIsAtRest( wrappedIsAtRest ) << 92 fWrappedProcessIsAtRest ( wrappedIsAtRest ),
70 fWrappedProcessIsAlong( wrappedIsAlongStep << 93 fWrappedProcessIsAlong ( wrappedIsAlongStep ),
71 fWrappedProcessIsPost( wrappedIsPostStep ) << 94 fWrappedProcessIsPost ( wrappedIsPostStep ),
>> 95 fWrappedProcessPostStepGPIL ( -1.0 ),
>> 96 fBiasingPostStepGPIL ( -1.0 ),
>> 97 fWrappedProcessInteractionLength ( -1.0 ),
>> 98 fWrappedProcessForceCondition ( NotForced ),
>> 99 fBiasingForceCondition ( NotForced ),
>> 100 fWrappedProcessAlongStepGPIL ( -1.0 ),
>> 101 fBiasingAlongStepGPIL ( -1.0 ),
>> 102 fWrappedProcessGPILSelection ( NotCandidateForSelection ),
>> 103 fBiasingGPILSelection ( NotCandidateForSelection ),
>> 104 fBiasingInteractionLaw ( nullptr ),
>> 105 fPreviousBiasingInteractionLaw ( nullptr ),
>> 106 fPhysicalInteractionLaw ( nullptr ),
>> 107 fOccurenceBiasingParticleChange ( nullptr ),
>> 108 fIamFirstGPIL ( false ),
>> 109 fProcessManager ( nullptr ),
>> 110 fSharedData ( nullptr )
72 { 111 {
73 for (G4int i = 0 ; i < 8 ; ++i) << 112 for (G4int i = 0 ; i < 8 ; i++) fFirstLastFlags[i] = false;
74 fFirstLastFlags[i] = false; <<
75 fResetInteractionLaws.Put( true ); 113 fResetInteractionLaws.Put( true );
76 fCommonStart.Put(true); 114 fCommonStart.Put(true);
77 fCommonEnd.Put(true); 115 fCommonEnd.Put(true);
78 fDoCommonConfigure.Put(true); 116 fDoCommonConfigure.Put(true);
79 117
80 SetProcessSubType(fWrappedProcess->GetProces 118 SetProcessSubType(fWrappedProcess->GetProcessSubType());
81 119
82 // -- create physical interaction law: 120 // -- create physical interaction law:
83 fPhysicalInteractionLaw = new G4InteractionL << 121 fPhysicalInteractionLaw = new G4InteractionLawPhysical("PhysicalInteractionLawFor("+GetProcessName()+")");
84 // -- instantiate particle change wrapper fo << 122 // -- instantiate particle change wrapper for occurence biaising:
85 fOccurenceBiasingParticleChange = new G4Part 123 fOccurenceBiasingParticleChange = new G4ParticleChangeForOccurenceBiasing("biasingPCfor"+GetProcessName());
86 // -- instantiate a "do nothing" particle ch 124 // -- instantiate a "do nothing" particle change:
87 fDummyParticleChange = new G4ParticleChangeF << 125 fDummyParticleChange = new G4ParticleChangeForNothing();
88 } 126 }
89 127
>> 128
>> 129
90 G4BiasingProcessInterface::~G4BiasingProcessIn 130 G4BiasingProcessInterface::~G4BiasingProcessInterface()
91 { 131 {
92 delete fPhysicalInteractionLaw; << 132 if ( fPhysicalInteractionLaw != 0 ) delete fPhysicalInteractionLaw;
93 delete fOccurenceBiasingParticleChange; << 133 if ( fOccurenceBiasingParticleChange ) delete fOccurenceBiasingParticleChange;
94 delete fDummyParticleChange; << 134 if ( fDummyParticleChange ) delete fDummyParticleChange;
95 } 135 }
96 136
97 const G4BiasingProcessSharedData* << 137
98 G4BiasingProcessInterface::GetSharedData( cons << 138 const G4BiasingProcessSharedData* G4BiasingProcessInterface::GetSharedData( const G4ProcessManager* mgr )
99 { 139 {
100 const auto & itr = G4BiasingProcessSharedDat << 140 G4MapCache< const G4ProcessManager*,
>> 141 G4BiasingProcessSharedData* >::const_iterator itr = G4BiasingProcessSharedData::fSharedDataMap.Find( mgr );
101 if ( itr != G4BiasingProcessSharedData::fSh 142 if ( itr != G4BiasingProcessSharedData::fSharedDataMap.End( ) )
102 { << 143 {
103 return (*itr).second; << 144 return (*itr).second;
104 } << 145 }
105 else return nullptr; << 146 else return 0;
106 } 147 }
107 148
108 149
109 void G4BiasingProcessInterface::StartTracking( 150 void G4BiasingProcessInterface::StartTracking(G4Track* track)
110 { 151 {
111 fCurrentTrack = track; 152 fCurrentTrack = track;
112 if ( fIsPhysicsBasedBiasing ) fWrappedProces 153 if ( fIsPhysicsBasedBiasing ) fWrappedProcess->StartTracking(fCurrentTrack);
113 fOccurenceBiasingOperation = nullpt << 154 fOccurenceBiasingOperation = 0;
114 fPreviousOccurenceBiasingOperation = nullpt << 155 fPreviousOccurenceBiasingOperation = 0;
115 fFinalStateBiasingOperation = nullpt << 156 fFinalStateBiasingOperation = 0;
116 fPreviousFinalStateBiasingOperation = nullpt << 157 fPreviousFinalStateBiasingOperation = 0;
117 fNonPhysicsBiasingOperation = nullpt << 158 fNonPhysicsBiasingOperation = 0;
118 fPreviousNonPhysicsBiasingOperation = nullpt << 159 fPreviousNonPhysicsBiasingOperation = 0;
119 fBiasingInteractionLaw = nullpt << 160 fBiasingInteractionLaw = 0;
120 fPreviousBiasingInteractionLaw = nullpt << 161 fPreviousBiasingInteractionLaw = 0;
121 162
122 fPreviousStepSize = -1.0; << 163 fPreviousStepSize = -1.0;
123 164
124 fResetWrappedProcessInteractionLength = fals 165 fResetWrappedProcessInteractionLength = false;
125 166
126 if ( fCommonStart.Get() ) 167 if ( fCommonStart.Get() )
127 { << 168 {
128 fCommonStart.Put( false );// = false; << 169 fCommonStart.Put( false );// = false;
129 fCommonEnd.Put( true );// = true; << 170 fCommonEnd.Put(true);// = true;
130 171
131 fSharedData->fCurrentBiasingOperator = nul << 172 fSharedData-> fCurrentBiasingOperator = 0;
132 fSharedData->fPreviousBiasingOperator = nu << 173 fSharedData->fPreviousBiasingOperator = 0;
133 174
134 // -- Add a "fSharedData->nStarting" here << 175 // -- §§ Add a "fSharedData->nStarting" here and outside bracket "fSharedData->nStarting++" and " if (fSharedData->nStarting) == fSharedData->(vector interface length)"
135 // -- call to the loop "StartTracking" of << 176 // -- §§ call to the loop "StartTracking" of operators"
136 177
137 for (std::size_t optr=0 ; optr<(G4VBiasing << 178 for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
138 { << 179 ( G4VBiasingOperator::GetBiasingOperators() )[optr]->StartTracking( fCurrentTrack );
139 (G4VBiasingOperator::GetBiasingOperators <<
140 } 180 }
141 } <<
142 } 181 }
143 182
>> 183
144 void G4BiasingProcessInterface::EndTracking() 184 void G4BiasingProcessInterface::EndTracking()
145 { 185 {
146 if ( fIsPhysicsBasedBiasing ) << 186 if ( fIsPhysicsBasedBiasing ) fWrappedProcess->EndTracking();
147 fWrappedProcess->EndTracking(); << 187 if ( fSharedData->fCurrentBiasingOperator) (fSharedData->fCurrentBiasingOperator)->ExitingBiasing( fCurrentTrack, this );
148 if ( fSharedData->fCurrentBiasingOperator) << 188 fBiasingInteractionLaw = 0;
149 (fSharedData->fCurrentBiasingOperator)->Ex <<
150 fBiasingInteractionLaw = nullptr; <<
151 189
152 // -- Inform operators of end of tracking: 190 // -- Inform operators of end of tracking:
153 if ( fCommonEnd.Get() ) 191 if ( fCommonEnd.Get() )
154 { <<
155 fCommonEnd .Put( false );// = false; <<
156 fCommonStart.Put( true );// = true; <<
157 <<
158 for ( std::size_t optr=0; optr<(G4VBiasing <<
159 { 192 {
160 (G4VBiasingOperator::GetBiasingOperators << 193 fCommonEnd.Put( false );// = false;
>> 194 fCommonStart.Put( true );// = true;
>> 195
>> 196 for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
>> 197 ( G4VBiasingOperator::GetBiasingOperators() )[optr]->EndTracking( );
>> 198
>> 199 // -- §§ for above loop, do as in StartTracking.
161 } 200 }
162 // -- for above loop, do as in StartTracki <<
163 } <<
164 } 201 }
165 202
166 G4double G4BiasingProcessInterface:: << 203
167 PostStepGetPhysicalInteractionLength( const G4 << 204
168 G4double << 205 G4double G4BiasingProcessInterface::PostStepGetPhysicalInteractionLength( const G4Track& track,
169 G4ForceC << 206 G4double previousStepSize,
>> 207 G4ForceCondition* condition )
170 { 208 {
>> 209
171 // ----------------------------------------- 210 // ---------------------------------------------------------------------------------------------------
172 // -- The "biasing process master" takes car 211 // -- The "biasing process master" takes care of updating the biasing operator, and for all biasing
173 // -- processes it invokes the PostStepGPIL 212 // -- processes it invokes the PostStepGPIL of physical wrapped processes (anticipate stepping manager
174 // -- call ! ) to make all cross-sections up 213 // -- call ! ) to make all cross-sections updated with current step, and hence available before the
175 // -- first call to the biasing operator. 214 // -- first call to the biasing operator.
176 // ----------------------------------------- 215 // ---------------------------------------------------------------------------------------------------
177 if ( fIamFirstGPIL ) 216 if ( fIamFirstGPIL )
178 { << 217 {
179 // -- Update previous biasing operator, an << 218 // -- Update previous biasing operator, and assume the operator stays the same by
180 // -- default and that it is not left at t << 219 // -- default and that it is not left at the beginning of this step. These
181 // -- assumptions might be wrong if there << 220 // -- assumptions might be wrong if there is a volume change (in paralllel or
182 // -- mass geometries) in what case the fl << 221 // -- mass geometries) in what case the flags will be updated.
183 fSharedData->fPreviousBiasingOperator = fS << 222 fSharedData->fPreviousBiasingOperator = fSharedData->fCurrentBiasingOperator;
184 fSharedData->fIsNewOperator = fa << 223 fSharedData->fIsNewOperator = false;
185 fSharedData->fLeavingPreviousOperator = fa << 224 fSharedData->fLeavingPreviousOperator = false;
186 // -- If new volume, either in mass or par << 225 // -- If new volume, either in mass or parallel geometries, get possible new biasing operator:
187 // --------------------------------------- << 226 // -------------------------------------------------------------------------------------------
188 // -- Get biasing operator in parallel geo << 227 // -- Get biasing operator in parallel geometries:
189 G4bool firstStepInParallelVolume = false; << 228 G4bool firstStepInParallelVolume = false;
190 if ( fSharedData->fParallelGeometriesLimit << 229 if ( fSharedData->fParallelGeometriesLimiterProcess )
191 { << 230 {
192 G4VBiasingOperator* newParallelOperator( << 231 G4VBiasingOperator* newParallelOperator( nullptr );
193 G4bool firstStep = ( track.GetCurrentSte << 232 G4bool firstStep = ( track.GetCurrentStepNumber() == 1 );
194 std::size_t iParallel = 0; << 233 size_t iParallel = 0;
195 for ( auto wasLimiting : fSharedData->fP << 234 for ( auto wasLimiting : fSharedData->fParallelGeometriesLimiterProcess->GetWasLimiting() )
196 { << 235 {
197 if ( firstStep || wasLimiting ) << 236 if ( firstStep || wasLimiting )
198 { << 237 {
199 firstStepInParallelVolume = true; << 238 firstStepInParallelVolume = true;
200 auto tmpParallelOperator = G4VBiasin << 239
201 GetBiasingOperator((fSharedData-> << 240 auto tmpParallelOperator = G4VBiasingOperator::GetBiasingOperator( (fSharedData->fParallelGeometriesLimiterProcess->GetCurrentVolumes()[iParallel])
202 ->GetCurrentVolumes()[iParallel]) << 241 ->GetLogicalVolume() );
203 ->GetLogicalVolume()); << 242 if ( newParallelOperator )
204 if ( newParallelOperator ) << 243 {
205 { << 244 if ( tmpParallelOperator )
206 if ( tmpParallelOperator ) << 245 {
207 { << 246 G4ExceptionDescription ed;
208 G4ExceptionDescription ed; << 247 ed << " Several biasing operators are defined at the same place in parallel geometries ! Found:\n";
209 ed << " Several biasing operator << 248 ed << " - `" << newParallelOperator->GetName() << "' and \n";
210 << " in parallel geometries ! << 249 ed << " - `" << tmpParallelOperator->GetName() << "'.\n";
211 ed << " - `" << newParallelOp << 250 ed << " Keeping `" << newParallelOperator->GetName() << "'. Behavior not guaranteed ! Please consider having only one operator at a place. " << G4endl;
212 ed << " - `" << tmpParallelOp << 251 G4Exception(" G4BiasingProcessInterface::PostStepGetPhysicalInteractionLength(...)",
213 ed << " Keeping `" << newParalle << 252 "BIAS.GEN.30",
214 << "'. Behavior not guarantee << 253 JustWarning,
215 << G4endl; << 254 ed);
216 G4Exception(" G4BiasingProcessIn << 255 }
217 "BIAS.GEN.30", JustW << 256 }
218 } << 257 else newParallelOperator = tmpParallelOperator;
219 } << 258 }
220 else newParallelOperator = tmpParall << 259 iParallel++;
221 } << 260 }
222 ++iParallel; << 261 fSharedData->fParallelGeometryOperator = newParallelOperator;
223 } << 262 } // -- end of " if ( fSharedData->fParallelGeometriesLimiterProcess )"
224 fSharedData->fParallelGeometryOperator = << 263
225 } // -- end of " if ( fSharedData->fParall << 264 // -- Get biasing operator in mass geometry:
226 << 265 // -- [§§ Note : bug with this first step ? Does not work if previous step was concurrently limited with geometry. Might make use of safety at last point ?]
227 // -- Get biasing operator in mass geometr << 266 G4bool firstStepInVolume = ( (track.GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) || (track.GetCurrentStepNumber() == 1) );
228 // -- [§§ Note : bug with this first ste << 267 // fSharedData->fIsNewOperator = false;
229 G4bool firstStepInVolume = ( (track.GetSt << 268 // fSharedData->fLeavingPreviousOperator = false;
230 || (track.GetCu << 269 if ( firstStepInVolume )
231 // fSharedData->fIsNewOperator << 270 {
232 // fSharedData->fLeavingPreviousOpera << 271 G4VBiasingOperator* newOperator = G4VBiasingOperator::GetBiasingOperator( track.GetVolume()->GetLogicalVolume() );
233 if ( firstStepInVolume ) << 272 fSharedData->fMassGeometryOperator = newOperator;
234 { << 273 if ( ( newOperator != nullptr ) && ( fSharedData->fParallelGeometryOperator != nullptr ) )
235 G4VBiasingOperator* newOperator = G4VBia << 274 {
236 GetBiasingOperator( track.GetVolume()- << 275 G4ExceptionDescription ed;
237 fSharedData->fMassGeometryOperator = new << 276 ed << " Biasing operators are defined at the same place in mass and parallel geometries ! Found:\n";
238 if ( ( newOperator != nullptr ) && ( fSh << 277 ed << " - `" << fSharedData->fParallelGeometryOperator->GetName() << "' in parallel geometry and \n";
239 { << 278 ed << " - `" << newOperator->GetName() << "' in mass geometry.\n";
240 G4ExceptionDescription ed; << 279 ed << " Keeping `" << fSharedData->fParallelGeometryOperator->GetName() << "'. Behavior not guaranteed ! Please consider having only one operator at a place. " << G4endl;
241 ed << " Biasing operators are defined << 280 G4Exception(" G4BiasingProcessInterface::PostStepGetPhysicalInteractionLength(...)",
242 ed << " - `" << fSharedData->fParal << 281 "BIAS.GEN.31",
243 ed << " - `" << newOperator->GetNam << 282 JustWarning,
244 ed << " Keeping `" << fSharedData->fPa << 283 ed);
245 G4Exception(" G4BiasingProcessInterfac << 284 }
246 "BIAS.GEN.31", JustWarning << 285 }
247 } << 286
248 } << 287 // -- conclude the operator selection, giving priority to parallel geometry (as told in exception message BIAS.GEN.30):
249 << 288 if ( firstStepInVolume || firstStepInParallelVolume )
250 // -- conclude the operator selection, giv << 289 {
251 if ( firstStepInVolume || firstStepInParal << 290 G4VBiasingOperator* newOperator = fSharedData->fParallelGeometryOperator;
252 { << 291 if ( newOperator == nullptr ) newOperator = fSharedData->fMassGeometryOperator;
253 G4VBiasingOperator* newOperator = fShare <<
254 if ( newOperator == nullptr ) <<
255 newOperator = fSharedData->fMassGeomet <<
256 292
257 fSharedData->fCurrentBiasingOperator = n << 293 fSharedData->fCurrentBiasingOperator = newOperator ;
258 294
259 if ( newOperator != fSharedData->fPrevio << 295 if ( newOperator != fSharedData->fPreviousBiasingOperator )
260 { << 296 {
261 fSharedData->fLeavingPreviousOperator << 297 fSharedData->fLeavingPreviousOperator = ( fSharedData->fPreviousBiasingOperator != nullptr ) ;
262 fSharedData->fIsNewOperator = ( newOpe << 298 fSharedData->fIsNewOperator = ( newOperator != nullptr );
263 } << 299 }
264 } << 300 }
265 << 301
266 // -- calls to wrapped process PostStepGPI << 302
267 // --------------------------------------- << 303 // -- calls to wrapped process PostStepGPIL's:
268 // -- Each physics wrapper process has its << 304 // -------------------------------------------
269 // -- fWrappedProcessPostStepGPIL , << 305 // -- Each physics wrapper process has its
270 // -- fWrappedProcessForceCondition , << 306 // -- fWrappedProcessPostStepGPIL ,
271 // -- fWrappedProcessInteractionLength << 307 // -- fWrappedProcessForceCondition ,
272 // -- updated. << 308 // -- fWrappedProcessInteractionLength
273 if ( fSharedData->fCurrentBiasingOperator << 309 // -- updated.
274 { << 310 if ( fSharedData->fCurrentBiasingOperator != nullptr )
275 for (std::size_t i=0; i<(fSharedData->fP << 311 {
276 { << 312 for ( size_t i = 0 ; i < (fSharedData->fPhysicsBiasingProcessInterfaces).size(); i++ )
277 (fSharedData->fPhysicsBiasingProcessIn << 313 (fSharedData->fPhysicsBiasingProcessInterfaces)[i]->InvokeWrappedProcessPostStepGPIL( track, previousStepSize, condition );
278 } << 314 }
279 } << 315 } // -- end of "if ( fIamFirstGPIL )"
280 } // -- end of "if ( fIamFirstGPIL )" << 316
>> 317
281 318
282 // -- Remember previous operator and propose 319 // -- Remember previous operator and proposed operations, if any, and reset:
283 // ----------------------------------------- 320 // -------------------------------------------------------------------------
284 // -- remember only in case some biasing mig 321 // -- remember only in case some biasing might be called
285 if ( ( fSharedData->fPreviousBiasingOperator << 322 if ( ( fSharedData->fPreviousBiasingOperator != 0 ) ||
286 ( fSharedData->fCurrentBiasingOperator << 323 ( fSharedData->fCurrentBiasingOperator != 0 ) )
287 { << 324 {
288 fPreviousOccurenceBiasingOperation = fOcc << 325 fPreviousOccurenceBiasingOperation = fOccurenceBiasingOperation;
289 fPreviousFinalStateBiasingOperation = fFin << 326 fPreviousFinalStateBiasingOperation = fFinalStateBiasingOperation;
290 fPreviousNonPhysicsBiasingOperation = fNon << 327 fPreviousNonPhysicsBiasingOperation = fNonPhysicsBiasingOperation;
291 fPreviousBiasingInteractionLaw = fBia << 328 fPreviousBiasingInteractionLaw = fBiasingInteractionLaw;
292 // -- reset: << 329 // -- reset:
293 fOccurenceBiasingOperation = null << 330 fOccurenceBiasingOperation = 0;
294 fFinalStateBiasingOperation = null << 331 fFinalStateBiasingOperation = 0;
295 fNonPhysicsBiasingOperation = null << 332 fNonPhysicsBiasingOperation = 0;
296 fBiasingInteractionLaw = null << 333 fBiasingInteractionLaw = 0;
297 // -- Physics PostStep and AlongStep GPIL << 334 // -- Physics PostStep and AlongStep GPIL
298 // fWrappedProcessPostStepGPIL : upda << 335 // fWrappedProcessPostStepGPIL : updated by InvokeWrappedProcessPostStepGPIL(...) above
299 fBiasingPostStepGPIL = DBL_ << 336 fBiasingPostStepGPIL = DBL_MAX;
300 // fWrappedProcessInteractionLength : upda << 337 // fWrappedProcessInteractionLength : updated by InvokeWrappedProcessPostStepGPIL(...) above; inverse of analog cross-section.
301 // fWrappedProcessForceCondition : upda << 338 // fWrappedProcessForceCondition : updated by InvokeWrappedProcessPostStepGPIL(...) above
302 fBiasingForceCondition = NotF << 339 fBiasingForceCondition = NotForced;
303 fWrappedProcessAlongStepGPIL = DBL_ << 340 fWrappedProcessAlongStepGPIL = DBL_MAX;
304 fBiasingAlongStepGPIL = DBL_ << 341 fBiasingAlongStepGPIL = DBL_MAX;
305 fWrappedProcessGPILSelection = NotC << 342 fWrappedProcessGPILSelection = NotCandidateForSelection;
306 fBiasingGPILSelection = NotC << 343 fBiasingGPILSelection = NotCandidateForSelection;
307 // -- for helper: << 344 // -- for helper:
308 fPreviousStepSize = prev << 345 fPreviousStepSize = previousStepSize;
309 } << 346 }
310 347
>> 348
311 // -- previous step size value; it is switch 349 // -- previous step size value; it is switched to zero if resetting a wrapped process:
312 // -- (same trick used than in InvokedWrappe 350 // -- (same trick used than in InvokedWrappedProcessPostStepGPIL )
313 G4double usedPreviousStepSize = previousStep 351 G4double usedPreviousStepSize = previousStepSize;
314 352
315 // ----------------------------------------- 353 // ----------------------------------------------
316 // -- If leaving a biasing operator, let it 354 // -- If leaving a biasing operator, let it know:
317 // ----------------------------------------- 355 // ----------------------------------------------
318 if ( fSharedData->fLeavingPreviousOperator ) 356 if ( fSharedData->fLeavingPreviousOperator )
319 { << 357 {
320 (fSharedData->fPreviousBiasingOperator)->E << 358 (fSharedData->fPreviousBiasingOperator)->ExitingBiasing( &track, this );
321 // -- if no further biasing operator, rese << 359 // -- if no further biasing operator, reset process behavior to standard tracking:
322 if ( fSharedData->fCurrentBiasingOperator << 360 if ( fSharedData->fCurrentBiasingOperator == 0 )
323 { << 361 {
324 ResetForUnbiasedTracking(); << 362 ResetForUnbiasedTracking();
325 if ( fIsPhysicsBasedBiasing ) << 363 if ( fIsPhysicsBasedBiasing )
326 { << 364 {
327 // -- if the physics process has been << 365 // -- if the physics process has been under occurence biasing, reset it:
328 if ( fResetWrappedProcessInteractionLe << 366 if ( fResetWrappedProcessInteractionLength )
329 { << 367 {
330 fResetWrappedProcessInteractionLengt << 368 fResetWrappedProcessInteractionLength = false;
331 fWrappedProcess->ResetNumberOfIntera << 369 fWrappedProcess->ResetNumberOfInteractionLengthLeft();
332 // -- We set "previous step size" as << 370 // -- We set "previous step size" as 0.0, to let the process believe this is first step:
333 usedPreviousStepSize = 0.0; << 371 usedPreviousStepSize = 0.0;
334 } << 372 }
335 } << 373 }
>> 374 }
336 } 375 }
337 } << 376
338 377
339 // ----------------------------------------- 378 // --------------------------------------------------------------
340 // -- no operator : analog tracking if physi 379 // -- no operator : analog tracking if physics-based, or nothing:
341 // ----------------------------------------- 380 // --------------------------------------------------------------
342 if ( fSharedData->fCurrentBiasingOperator == << 381 if ( fSharedData->fCurrentBiasingOperator == 0 )
343 { <<
344 // -- take note of the "usedPreviousStepSi <<
345 if ( fIsPhysicsBasedBiasing ) <<
346 { 382 {
347 return fWrappedProcess->PostStepGetPhysi << 383 // -- take note of the "usedPreviousStepSize" value:
>> 384 if ( fIsPhysicsBasedBiasing ) return fWrappedProcess->PostStepGetPhysicalInteractionLength(track, usedPreviousStepSize, condition);
>> 385 else
>> 386 {
>> 387 *condition = NotForced;
>> 388 return DBL_MAX;
>> 389 }
348 } 390 }
349 else << 391
350 { << 392
351 *condition = NotForced; <<
352 return DBL_MAX; <<
353 } <<
354 } <<
355 393
356 // ----------------------------------------- 394 // --------------------------------------------------
357 // -- A biasing operator exists. Proceed wit 395 // -- A biasing operator exists. Proceed with
358 // -- treating non-physics and physics biasi 396 // -- treating non-physics and physics biasing cases:
359 //------------------------------------------ 397 //---------------------------------------------------
360 398
361 // -- non-physics-based biasing case: 399 // -- non-physics-based biasing case:
362 // ---------------------------------- 400 // ----------------------------------
363 if ( !fIsPhysicsBasedBiasing ) 401 if ( !fIsPhysicsBasedBiasing )
364 { << 402 {
365 fNonPhysicsBiasingOperation = (fSharedData << 403 fNonPhysicsBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedNonPhysicsBiasingOperation( &track, this );
366 if ( fNonPhysicsBiasingOperation == nullpt << 404 if ( fNonPhysicsBiasingOperation == 0 )
367 { << 405 {
368 *condition = NotForced; << 406 *condition = NotForced;
369 return DBL_MAX; << 407 return DBL_MAX;
>> 408 }
>> 409 return fNonPhysicsBiasingOperation->DistanceToApplyOperation(&track, previousStepSize, condition);
370 } 410 }
371 return fNonPhysicsBiasingOperation->Distan << 411
372 } <<
373 412
374 // -- Physics based biasing case: 413 // -- Physics based biasing case:
375 // ------------------------------ 414 // ------------------------------
376 // -- Ask for possible GPIL biasing operatio 415 // -- Ask for possible GPIL biasing operation:
377 fOccurenceBiasingOperation = (fSharedData->f 416 fOccurenceBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedOccurenceBiasingOperation( &track, this );
378 417
379 // -- no operation for occurrence biasing, a << 418
380 if ( fOccurenceBiasingOperation == nullptr ) << 419 // -- no operation for occurence biasing, analog GPIL returns the wrapped process GPIL and condition values
381 { << 420 if ( fOccurenceBiasingOperation == 0 )
382 *condition = fWrappedProcessForceCondition << 421 {
383 return fWrappedProcessPostStepGPIL; << 422 *condition = fWrappedProcessForceCondition;
384 } << 423 return fWrappedProcessPostStepGPIL;
>> 424 }
385 425
386 // -- A valid GPIL biasing operation has bee 426 // -- A valid GPIL biasing operation has been proposed:
387 // -- 0) remember wrapped process will need 427 // -- 0) remember wrapped process will need to be reset on biasing exit, if particle survives:
388 fResetWrappedProcessInteractionLength = true 428 fResetWrappedProcessInteractionLength = true;
389 // -- 1) update process interaction length f 429 // -- 1) update process interaction length for reference analog interaction law ( fWrappedProcessInteractionLength updated/collected above):
390 fPhysicalInteractionLaw->SetPhysicalCrossSec 430 fPhysicalInteractionLaw->SetPhysicalCrossSection( 1.0 / fWrappedProcessInteractionLength );
391 // -- 2) Collect biasing interaction law: 431 // -- 2) Collect biasing interaction law:
392 // -- The interaction law pointer is coll 432 // -- The interaction law pointer is collected as a const pointer to the interaction law object.
393 // -- This interaction law will be kept u 433 // -- This interaction law will be kept under control of the biasing operation, which is the only
394 // -- entity that will change the state o 434 // -- entity that will change the state of the biasing interaction law.
395 // -- The force condition for biasing is 435 // -- The force condition for biasing is asked at the same time, passing the analog one as default:
396 fBiasingForceCondition = fWrappedProcessForc << 436 fBiasingForceCondition = fWrappedProcessForceCondition;
397 fBiasingInteractionLaw = fOccurenceBiasingOp << 437 fBiasingInteractionLaw = fOccurenceBiasingOperation->ProvideOccurenceBiasingInteractionLaw( this, fBiasingForceCondition );
398 // -- 3) Ask operation to sample the biasing 438 // -- 3) Ask operation to sample the biasing interaction law:
399 fBiasingPostStepGPIL = fBiasingInteractionLa << 439 fBiasingPostStepGPIL = fBiasingInteractionLaw->GetSampledInteractionLength();
400 440
401 // -- finish 441 // -- finish
402 *condition = fBiasingForceCondition; 442 *condition = fBiasingForceCondition;
403 return fBiasingPostStepGPIL; << 443 return fBiasingPostStepGPIL;
>> 444
404 } 445 }
405 446
>> 447
>> 448
406 G4VParticleChange* G4BiasingProcessInterface:: 449 G4VParticleChange* G4BiasingProcessInterface::PostStepDoIt(const G4Track& track,
407 << 450 const G4Step& step)
408 { 451 {
>> 452
409 // --------------------------------------- 453 // ---------------------------------------
410 // -- case outside of volume with biasing: 454 // -- case outside of volume with biasing:
411 // --------------------------------------- 455 // ---------------------------------------
412 if ( fSharedData->fCurrentBiasingOperator == << 456 if ( fSharedData->fCurrentBiasingOperator == 0 ) return fWrappedProcess->PostStepDoIt(track, step);
413 return fWrappedProcess->PostStepDoIt(track <<
414 457
415 // ---------------------------- 458 // ----------------------------
416 // -- non-physics biasing case: 459 // -- non-physics biasing case:
417 // ---------------------------- 460 // ----------------------------
418 if ( !fIsPhysicsBasedBiasing ) 461 if ( !fIsPhysicsBasedBiasing )
419 { << 462 {
420 G4VParticleChange* particleChange = fNonPh << 463 G4VParticleChange* particleChange = fNonPhysicsBiasingOperation->GenerateBiasingFinalState( &track, &step );
421 (fSharedData->fCurrentBiasingOperator)->Re << 464 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC_NonPhysics, fNonPhysicsBiasingOperation, particleChange );
422 return particleChange; << 465 return particleChange;
423 } << 466 }
424 467
425 // -- physics biasing case: 468 // -- physics biasing case:
426 // ------------------------ 469 // ------------------------
427 // -- It proceeds with the following logic: 470 // -- It proceeds with the following logic:
428 // -- 1) Obtain the final state 471 // -- 1) Obtain the final state
429 // -- This final state may be analog or b 472 // -- This final state may be analog or biased.
430 // -- The biased final state is obtained 473 // -- The biased final state is obtained through a biasing operator
431 // -- returned by the operator. 474 // -- returned by the operator.
432 // -- 2) The biased final state may be asked 475 // -- 2) The biased final state may be asked to be "force as it is"
433 // -- in what case the particle change is 476 // -- in what case the particle change is returned as is to the
434 // -- stepping. 477 // -- stepping.
435 // -- In all other cases (analog final st 478 // -- In all other cases (analog final state or biased final but
436 // -- not forced) the final state weight 479 // -- not forced) the final state weight may be modified by the
437 // -- occurrence biasing, if such an occu << 480 // -- occurence biasing, if such an occurence biasing is at play.
438 481
439 // -- Get final state, biased or analog: 482 // -- Get final state, biased or analog:
440 G4VParticleChange* finalStateParticleChange; << 483 G4VParticleChange* finalStateParticleChange;
441 G4BiasingAppliedCase BAC; 484 G4BiasingAppliedCase BAC;
442 fFinalStateBiasingOperation = (fSharedData-> 485 fFinalStateBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedFinalStateBiasingOperation( &track, this );
443 // -- Flag below is to force the biased gene << 486 // -- Flag below is to force the biased generated particle change to be retruned "as is" to the stepping, disregarding there
444 // -- was or not a occurrence biasing that w << 487 // -- was or not a occurence biasing that would apply. Weight relevance under full responsibility of the biasing operation.
445 G4bool forceBiasedFinalState = false; 488 G4bool forceBiasedFinalState = false;
446 if ( fFinalStateBiasingOperation != nullptr << 489 if ( fFinalStateBiasingOperation != 0 )
447 { << 490 {
448 finalStateParticleChange = fFinalStateBias << 491 finalStateParticleChange = fFinalStateBiasingOperation->ApplyFinalStateBiasing( this, &track, &step, forceBiasedFinalState );
449 BAC = BAC_FinalState; << 492 BAC = BAC_FinalState;
450 } << 493 }
451 else 494 else
452 { << 495 {
453 finalStateParticleChange = fWrappedProcess << 496 finalStateParticleChange = fWrappedProcess->PostStepDoIt(track, step);
454 BAC = BAC_None ; << 497 BAC = BAC_None ;
455 } << 498 }
456 << 499
457 // -- if no occurrence biasing operation, we << 500 // -- if no occurence biasing operation, we're done:
458 if ( fOccurenceBiasingOperation == nullptr ) << 501 if ( fOccurenceBiasingOperation == 0 )
459 { << 502 {
460 (fSharedData->fCurrentBiasingOperator)->Re << 503 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC, fFinalStateBiasingOperation, finalStateParticleChange );
461 return finalStateParticleChange; << 504 return finalStateParticleChange;
462 } << 505 }
463 506
464 // -- if biased final state has been asked t 507 // -- if biased final state has been asked to be forced, we're done:
465 if ( forceBiasedFinalState ) 508 if ( forceBiasedFinalState )
466 { << 509 {
467 (fSharedData->fCurrentBiasingOperator)->Re << 510 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC, fFinalStateBiasingOperation, finalStateParticleChange );
468 return finalStateParticleChange; << 511 return finalStateParticleChange;
469 } << 512 }
>> 513
470 514
471 // -- If occurrence biasing, applies the occ << 515 // -- If occurence biasing, applies the occurence biasing weight correction on top of final state (biased or not):
472 G4double weightForInteraction = 1.0; 516 G4double weightForInteraction = 1.0;
473 if ( !fBiasingInteractionLaw->IsSingular() ) << 517 if ( !fBiasingInteractionLaw->IsSingular() ) weightForInteraction =
474 { << 518 fPhysicalInteractionLaw->ComputeEffectiveCrossSectionAt(step.GetStepLength()) /
475 weightForInteraction = fPhysicalInteractio << 519 fBiasingInteractionLaw ->ComputeEffectiveCrossSectionAt(step.GetStepLength());
476 / fBiasingInteraction <<
477 } <<
478 else 520 else
479 { << 521 {
480 // -- at this point effective XS can only << 522 // -- at this point effective XS can only be infinite, if not, there is a logic problem
481 if ( !fBiasingInteractionLaw->IsEffectiveC << 523 if ( !fBiasingInteractionLaw->IsEffectiveCrossSectionInfinite() )
>> 524 {
>> 525 G4ExceptionDescription ed;
>> 526 ed << "Internal inconsistency in cross-section handling. Please report !" << G4endl;
>> 527 G4Exception(" G4BiasingProcessInterface::PostStepDoIt(...)",
>> 528 "BIAS.GEN.02",
>> 529 JustWarning,
>> 530 ed);
>> 531 // -- if XS is infinite, weight is zero (and will stay zero), but we'll do differently.
>> 532 // -- Should foresee in addition something to remember that in case of singular
>> 533 // -- distribution, weight can only be partly calculated
>> 534 }
>> 535 }
>> 536
>> 537 if ( weightForInteraction <= 0. )
482 { 538 {
483 G4ExceptionDescription ed; 539 G4ExceptionDescription ed;
484 ed << "Internal inconsistency in cross-s << 540 ed << " Negative interaction weight : w_I = "
>> 541 << weightForInteraction <<
>> 542 " XS_I(phys) = " << fBiasingInteractionLaw ->ComputeEffectiveCrossSectionAt(step.GetStepLength()) <<
>> 543 " XS_I(bias) = " << fPhysicalInteractionLaw->ComputeEffectiveCrossSectionAt(step.GetStepLength()) <<
>> 544 " step length = " << step.GetStepLength() <<
>> 545 " Interaction law = `" << fBiasingInteractionLaw << "'" <<
>> 546 G4endl;
485 G4Exception(" G4BiasingProcessInterface: 547 G4Exception(" G4BiasingProcessInterface::PostStepDoIt(...)",
486 "BIAS.GEN.02", JustWarning, << 548 "BIAS.GEN.03",
487 // -- if XS is infinite, weight is zero << 549 JustWarning,
488 // -- Should foresee in addition somethi << 550 ed);
489 // -- distribution, weight can only be p << 551
490 } 552 }
491 } <<
492 553
493 if ( weightForInteraction <= 0. ) << 554 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC,
494 { << 555 fOccurenceBiasingOperation, weightForInteraction,
495 G4ExceptionDescription ed; << 556 fFinalStateBiasingOperation, finalStateParticleChange );
496 ed << " Negative interaction weight : w_I << 557
497 << weightForInteraction << " XS_I(phys <<
498 << fBiasingInteractionLaw ->ComputeEffe <<
499 <<" XS_I(bias) = " <<
500 << fPhysicalInteractionLaw->ComputeEffe <<
501 << " step length = " << step.GetStepLen <<
502 << " Interaction law = `" << fBiasingIn <<
503 << G4endl; <<
504 G4Exception(" G4BiasingProcessInterface::P <<
505 "BIAS.GEN.03", JustWarning, ed <<
506 } <<
507 <<
508 (fSharedData->fCurrentBiasingOperator) <<
509 ->ReportOperationApplied( this, BAC, fOcc <<
510 weightForIntera <<
511 fFinalStateBias <<
512 finalStateParti <<
513 558
514 fOccurenceBiasingParticleChange->SetOccurenc 559 fOccurenceBiasingParticleChange->SetOccurenceWeightForInteraction( weightForInteraction );
515 fOccurenceBiasingParticleChange->SetSecondar 560 fOccurenceBiasingParticleChange->SetSecondaryWeightByProcess( true );
516 fOccurenceBiasingParticleChange->SetWrappedP 561 fOccurenceBiasingParticleChange->SetWrappedParticleChange( finalStateParticleChange );
517 fOccurenceBiasingParticleChange->ProposeTrac 562 fOccurenceBiasingParticleChange->ProposeTrackStatus( finalStateParticleChange->GetTrackStatus() );
518 fOccurenceBiasingParticleChange->StealSecond 563 fOccurenceBiasingParticleChange->StealSecondaries(); // -- this also makes weightForInteraction applied to secondaries stolen
519 564
520 // -- finish: 565 // -- finish:
521 return fOccurenceBiasingParticleChange; 566 return fOccurenceBiasingParticleChange;
>> 567
522 } 568 }
523 569
>> 570
524 // -- AlongStep methods: 571 // -- AlongStep methods:
525 G4double G4BiasingProcessInterface:: << 572 G4double G4BiasingProcessInterface::AlongStepGetPhysicalInteractionLength(const G4Track& track,
526 AlongStepGetPhysicalInteractionLength(const G4 << 573 G4double previousStepSize,
527 G4 << 574 G4double currentMinimumStep,
528 G4 << 575 G4double& proposedSafety,
529 G4 << 576 G4GPILSelection* selection)
530 G4 <<
531 { 577 {
>> 578
532 // -- for helper methods: 579 // -- for helper methods:
533 fCurrentMinimumStep = currentMinimumStep; 580 fCurrentMinimumStep = currentMinimumStep;
534 fProposedSafety = proposedSafety; << 581 fProposedSafety = proposedSafety;
>> 582
535 583
536 // -- initialization default case: 584 // -- initialization default case:
537 fWrappedProcessAlongStepGPIL = DBL_MAX; << 585 fWrappedProcessAlongStepGPIL = DBL_MAX;
538 *selection = NotCandidateF 586 *selection = NotCandidateForSelection;
539 // --------------------------------------- 587 // ---------------------------------------
540 // -- case outside of volume with biasing: 588 // -- case outside of volume with biasing:
541 // --------------------------------------- 589 // ---------------------------------------
542 if ( fSharedData->fCurrentBiasingOperator == << 590 if ( fSharedData->fCurrentBiasingOperator == 0 )
543 { << 591 {
544 if ( fWrappedProcessIsAlong ) << 592 if ( fWrappedProcessIsAlong ) fWrappedProcessAlongStepGPIL =
545 fWrappedProcessAlongStepGPIL = fWrappedP << 593 fWrappedProcess->AlongStepGetPhysicalInteractionLength(track,
546 ->AlongStepGetPhysicalInteractionLengt << 594 previousStepSize,
547 << 595 currentMinimumStep,
548 << 596 proposedSafety,
549 return fWrappedProcessAlongStepGPIL; << 597 selection);
550 } << 598 return fWrappedProcessAlongStepGPIL;
>> 599 }
551 600
552 // ----------------------------------------- 601 // --------------------------------------------------------------------
553 // -- non-physics based biasing: no along op 602 // -- non-physics based biasing: no along operation expected (for now):
554 // ----------------------------------------- 603 // --------------------------------------------------------------------
555 if ( !fIsPhysicsBasedBiasing ) return fWrapp 604 if ( !fIsPhysicsBasedBiasing ) return fWrappedProcessAlongStepGPIL;
556 605
557 // ---------------------- 606 // ----------------------
558 // -- physics-based case: 607 // -- physics-based case:
559 // ---------------------- 608 // ----------------------
560 if ( fOccurenceBiasingOperation == nullptr ) << 609 if ( fOccurenceBiasingOperation == 0 )
561 { << 610 {
562 if ( fWrappedProcessIsAlong ) << 611 if ( fWrappedProcessIsAlong ) fWrappedProcessAlongStepGPIL =
563 fWrappedProcessAlongStepGPIL = fWrappedP << 612 fWrappedProcess->AlongStepGetPhysicalInteractionLength(track,
564 ->AlongStepGetPhysicalInteractionLengt << 613 previousStepSize,
565 << 614 currentMinimumStep,
566 << 615 proposedSafety,
567 return fWrappedProcessAlongStepGPIL; << 616 selection);
568 } << 617 return fWrappedProcessAlongStepGPIL;
569 << 618 }
570 // ----------------------------------------- << 619
571 // -- From here we have an valid occurrence << 620
572 // ----------------------------------------- << 621 // ----------------------------------------------------------
>> 622 // -- From here we have an valid occurence biasing operation:
>> 623 // ----------------------------------------------------------
573 // -- Give operation opportunity to shorten 624 // -- Give operation opportunity to shorten step proposed by physics process:
574 fBiasingAlongStepGPIL = fOccurenceBiasingOpe << 625 fBiasingAlongStepGPIL = fOccurenceBiasingOperation->ProposeAlongStepLimit( this );
575 G4double minimumStep = fBiasingAlongStepGPI << 626 G4double minimumStep = fBiasingAlongStepGPIL < currentMinimumStep ? fBiasingAlongStepGPIL : currentMinimumStep ;
576 ? fBiasingAlongStepGPI <<
577 // -- wrapped process is called with minimum 627 // -- wrapped process is called with minimum step ( <= currentMinimumStep passed ) : an along process can not
578 // -- have its operation stretched over what 628 // -- have its operation stretched over what it expects:
579 if ( fWrappedProcessIsAlong ) 629 if ( fWrappedProcessIsAlong )
580 { << 630 {
581 fWrappedProcessAlongStepGPIL = fWrappedPro << 631 fWrappedProcessAlongStepGPIL = fWrappedProcess->AlongStepGetPhysicalInteractionLength(track,
582 ->AlongStepGetPhysicalInteractionLength( << 632 previousStepSize,
583 << 633 minimumStep,
584 << 634 proposedSafety,
585 fWrappedProcessGPILSelection = *selection; << 635 selection);
586 fBiasingGPILSelection = fOccurenceBiasingO << 636 fWrappedProcessGPILSelection = *selection;
587 ->ProposeGPILSelection( fWrappedProcessG << 637 fBiasingGPILSelection = fOccurenceBiasingOperation->ProposeGPILSelection( fWrappedProcessGPILSelection );
588 } << 638 }
589 else 639 else
590 { << 640 {
591 fBiasingGPILSelection = fOccurenceBiasingO << 641 fBiasingGPILSelection = fOccurenceBiasingOperation->ProposeGPILSelection( NotCandidateForSelection );
592 ->ProposeGPILSelection( NotCandidateForS << 642 fWrappedProcessAlongStepGPIL = fBiasingAlongStepGPIL;
593 fWrappedProcessAlongStepGPIL = fBiasingAlo << 643 }
594 } <<
595 644
596 *selection = fBiasingGPILSelection; 645 *selection = fBiasingGPILSelection;
597 646
598 return fWrappedProcessAlongStepGPIL; 647 return fWrappedProcessAlongStepGPIL;
>> 648
599 } 649 }
600 650
601 G4VParticleChange* << 651 G4VParticleChange* G4BiasingProcessInterface::AlongStepDoIt(const G4Track& track,
602 G4BiasingProcessInterface::AlongStepDoIt(const << 652 const G4Step& step)
603 const <<
604 { 653 {
>> 654
605 // --------------------------------------- 655 // ---------------------------------------
606 // -- case outside of volume with biasing: 656 // -- case outside of volume with biasing:
607 // --------------------------------------- 657 // ---------------------------------------
608 if ( fSharedData->fCurrentBiasingOperator == << 658 if ( fSharedData->fCurrentBiasingOperator == 0 )
609 { <<
610 if ( fWrappedProcessIsAlong ) <<
611 { 659 {
612 return fWrappedProcess->AlongStepDoIt(tr << 660 if ( fWrappedProcessIsAlong ) return fWrappedProcess->AlongStepDoIt(track, step);
>> 661 else
>> 662 {
>> 663 fDummyParticleChange->Initialize( track );
>> 664 return fDummyParticleChange;
>> 665 }
613 } 666 }
614 else <<
615 { <<
616 fDummyParticleChange->Initialize( track <<
617 return fDummyParticleChange; <<
618 } <<
619 } <<
620 667
621 // ----------------------------------- 668 // -----------------------------------
622 // -- case inside volume with biasing: 669 // -- case inside volume with biasing:
623 // ----------------------------------- 670 // -----------------------------------
624 if ( fWrappedProcessIsAlong ) << 671 if ( fWrappedProcessIsAlong ) fOccurenceBiasingParticleChange->SetWrappedParticleChange( fWrappedProcess->AlongStepDoIt(track, step) );
625 { <<
626 fOccurenceBiasingParticleChange <<
627 ->SetWrappedParticleChange(fWrappedProce <<
628 } <<
629 else 672 else
630 { << 673 {
631 fOccurenceBiasingParticleChange->SetWrappe << 674 fOccurenceBiasingParticleChange->SetWrappedParticleChange ( 0 );
632 fOccurenceBiasingParticleChange->ProposeTr << 675 fOccurenceBiasingParticleChange->ProposeTrackStatus( track.GetTrackStatus() );
633 } << 676 }
634 G4double weightForNonInteraction (1.0); 677 G4double weightForNonInteraction (1.0);
635 if ( fBiasingInteractionLaw != nullptr ) << 678 if ( fBiasingInteractionLaw != 0 )
636 { << 679 {
637 weightForNonInteraction = << 680 weightForNonInteraction =
638 fPhysicalInteractionLaw->ComputeNonInter << 681 fPhysicalInteractionLaw->ComputeNonInteractionProbabilityAt(step.GetStepLength()) /
639 fBiasingInteractionLaw ->ComputeNonInter << 682 fBiasingInteractionLaw ->ComputeNonInteractionProbabilityAt(step.GetStepLength());
640 683
641 fOccurenceBiasingOperation->AlongMoveBy( t << 684 fOccurenceBiasingOperation->AlongMoveBy( this, &step, weightForNonInteraction );
642 685
643 if ( weightForNonInteraction <= 0. ) << 686 if ( weightForNonInteraction <= 0. )
644 { << 687 {
645 G4ExceptionDescription ed; << 688 G4ExceptionDescription ed;
646 ed << " Negative non interaction weight << 689 ed << " Negative non interaction weight : w_NI = " << weightForNonInteraction <<
647 " p_NI(phys) = " << fPhysicalInte << 690 " p_NI(phys) = " << fPhysicalInteractionLaw->ComputeNonInteractionProbabilityAt(step.GetStepLength()) <<
648 " p_NI(bias) = " << fBiasingInter << 691 " p_NI(bias) = " << fBiasingInteractionLaw ->ComputeNonInteractionProbabilityAt(step.GetStepLength()) <<
649 " step length = " << step.GetSte << 692 " step length = " << step.GetStepLength() <<
650 " biasing interaction law = `" << << 693 " biasing interaction law = `" << fBiasingInteractionLaw->GetName() << "'" << G4endl;
651 G4Exception(" G4BiasingProcessInterface: << 694 G4Exception(" G4BiasingProcessInterface::AlongStepDoIt(...)",
652 "BIAS.GEN.04", JustWarning, << 695 "BIAS.GEN.04",
>> 696 JustWarning,
>> 697 ed);
>> 698 }
>> 699
653 } 700 }
654 } <<
655 701
656 fOccurenceBiasingParticleChange << 702 fOccurenceBiasingParticleChange->SetOccurenceWeightForNonInteraction( weightForNonInteraction );
657 ->SetOccurenceWeightForNonInteraction( wei <<
658 703
659 return fOccurenceBiasingParticleChange; 704 return fOccurenceBiasingParticleChange;
>> 705
660 } 706 }
661 707
662 // -- AtRest methods 708 // -- AtRest methods
663 G4double G4BiasingProcessInterface:: << 709 G4double G4BiasingProcessInterface::AtRestGetPhysicalInteractionLength(const G4Track& track,
664 AtRestGetPhysicalInteractionLength(const G4Tra << 710 G4ForceCondition* condition)
665 G4ForceCond <<
666 { 711 {
667 return fWrappedProcess->AtRestGetPhysicalInt << 712 return fWrappedProcess->AtRestGetPhysicalInteractionLength(track, condition);
668 } 713 }
669 <<
670 G4VParticleChange* G4BiasingProcessInterface:: 714 G4VParticleChange* G4BiasingProcessInterface::AtRestDoIt(const G4Track& track,
671 << 715 const G4Step& step)
672 { 716 {
673 return fWrappedProcess->AtRestDoIt(track, st << 717 return fWrappedProcess->AtRestDoIt(track, step);
674 } 718 }
675 719
676 G4bool G4BiasingProcessInterface::IsApplicable << 720
>> 721 G4bool G4BiasingProcessInterface::IsApplicable(const G4ParticleDefinition& pd)
677 { 722 {
678 if ( fWrappedProcess != nullptr ) << 723 if ( fWrappedProcess != 0 ) return fWrappedProcess->IsApplicable(pd);
679 return fWrappedProcess->IsApplicable(pd); << 724 else return true;
680 else <<
681 return true; <<
682 } 725 }
683 726
684 void G4BiasingProcessInterface::SetMasterProc << 727
>> 728 void G4BiasingProcessInterface::SetMasterProcess(G4VProcess* masterP)
685 { 729 {
686 // -- Master for this process: 730 // -- Master for this process:
687 G4VProcess::SetMasterProcess(masterP); 731 G4VProcess::SetMasterProcess(masterP);
688 // -- Master for wrapped process: 732 // -- Master for wrapped process:
689 if ( fWrappedProcess != nullptr ) << 733 if ( fWrappedProcess != 0 )
690 { << 734 {
691 const G4BiasingProcessInterface* thisWrapp << 735 const G4BiasingProcessInterface* thisWrapperMaster = (const G4BiasingProcessInterface *)GetMasterProcess();
692 = (const G4BiasingProcessInterface *)Get << 736 // -- paranoia check: (?)
693 // -- paranoia check: (?) << 737 G4VProcess* wrappedMaster = 0;
694 G4VProcess* wrappedMaster = nullptr; << 738 wrappedMaster = thisWrapperMaster->GetWrappedProcess();
695 wrappedMaster = thisWrapperMaster->GetWrap << 739 fWrappedProcess->SetMasterProcess( wrappedMaster );
696 fWrappedProcess->SetMasterProcess( wrapped << 740 }
697 } <<
698 } 741 }
699 742
700 void G4BiasingProcessInterface:: << 743
701 BuildPhysicsTable(const G4ParticleDefinition& << 744 void G4BiasingProcessInterface::BuildPhysicsTable(const G4ParticleDefinition& pd)
702 { 745 {
703 // -- Sequential mode : called second (after 746 // -- Sequential mode : called second (after PreparePhysicsTable(..))
704 // -- MT mode : called second (after 747 // -- MT mode : called second (after PreparePhysicsTable(..)) by master thread.
705 // -- Corresponding proces 748 // -- Corresponding process instance not used then by tracking.
706 // -- PreparePhysicsTable(...) has been call 749 // -- PreparePhysicsTable(...) has been called first for all processes,
707 // -- so the first/last flags and G4BiasingP 750 // -- so the first/last flags and G4BiasingProcessInterface vector of processes have
708 // -- been properly setup, fIamFirstGPIL is 751 // -- been properly setup, fIamFirstGPIL is valid.
709 if ( fWrappedProcess != nullptr ) << 752 if ( fWrappedProcess != 0 )
710 { << 753 {
711 fWrappedProcess->BuildPhysicsTable(pd); << 754 fWrappedProcess->BuildPhysicsTable(pd);
712 } << 755 }
713 756
714 if ( fIamFirstGPIL ) 757 if ( fIamFirstGPIL )
715 { << 758 {
716 // -- Re-order vector of processes to matc << 759 // -- Re-order vector of processes to match that of the GPIL
717 // -- (made for fIamFirstGPIL, but importa << 760 // -- (made for fIamFirstGPIL, but important is to have it made once):
718 ReorderBiasingVectorAsGPIL(); << 761 ReorderBiasingVectorAsGPIL();
719 // -- Let operators to configure themselve << 762 // -- Let operators to configure themselves for the master thread or for sequential mode.
720 // -- Intended here is in particular the r << 763 // -- Intended here is in particular the registration to physics model catalog.
721 // -- The fDoCommonConfigure is to ensure << 764 // -- The fDoCommonConfigure is to ensure that this Configure is made by only one process (othewise each first process makes the call):
722 if ( fDoCommonConfigure.Get() ) << 765 if ( fDoCommonConfigure.Get() )
723 { << 766 {
724 for ( std::size_t optr=0; optr<(G4VBiasi << 767 for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
725 { << 768 ( G4VBiasingOperator::GetBiasingOperators() )[optr]->Configure( );
726 (G4VBiasingOperator::GetBiasingOperato << 769 fDoCommonConfigure.Put(false);
727 } << 770 }
728 fDoCommonConfigure.Put(false); << 771
729 } 772 }
730 } <<
731 } 773 }
732 774
733 void G4BiasingProcessInterface:: << 775
734 PreparePhysicsTable(const G4ParticleDefinition << 776 void G4BiasingProcessInterface::PreparePhysicsTable(const G4ParticleDefinition& pd)
735 { 777 {
736 // -- Sequential mode : called first (before 778 // -- Sequential mode : called first (before BuildPhysicsTable(..))
737 // -- MT mode : called first (before 779 // -- MT mode : called first (before BuildPhysicsTable(..)) by master thread.
738 // -- Corresponding proces 780 // -- Corresponding process instance not used then by tracking.
739 // -- Let process finding its first/last pos 781 // -- Let process finding its first/last position in the process manager:
740 SetUpFirstLastFlags(); 782 SetUpFirstLastFlags();
741 if ( fWrappedProcess != nullptr ) << 783 if ( fWrappedProcess != 0 )
742 { << 784 {
743 fWrappedProcess->PreparePhysicsTable(pd); << 785 fWrappedProcess->PreparePhysicsTable(pd);
744 } << 786 }
745 } 787 }
746 788
747 G4bool G4BiasingProcessInterface:: << 789
748 StorePhysicsTable(const G4ParticleDefinition* << 790 G4bool G4BiasingProcessInterface::StorePhysicsTable(const G4ParticleDefinition* pd, const G4String& s, G4bool f)
749 { 791 {
750 if ( fWrappedProcess != nullptr ) << 792 if ( fWrappedProcess != 0 ) return fWrappedProcess->StorePhysicsTable(pd, s, f);
751 return fWrappedProcess->StorePhysicsTable( << 793 else return false;
752 else <<
753 return false; <<
754 } 794 }
755 795
756 G4bool G4BiasingProcessInterface:: << 796
757 RetrievePhysicsTable(const G4ParticleDefinitio << 797 G4bool G4BiasingProcessInterface::RetrievePhysicsTable(const G4ParticleDefinition* pd, const G4String& s, G4bool f)
758 { 798 {
759 if ( fWrappedProcess != nullptr ) << 799 if ( fWrappedProcess != 0 ) return fWrappedProcess->RetrievePhysicsTable(pd, s, f);
760 return fWrappedProcess->RetrievePhysicsTab << 800 else return false;
761 else <<
762 return false; <<
763 } 801 }
764 802
>> 803
765 void G4BiasingProcessInterface::SetProcessMana 804 void G4BiasingProcessInterface::SetProcessManager(const G4ProcessManager* mgr)
766 { 805 {
767 if ( fWrappedProcess != nullptr ) << 806 if ( fWrappedProcess != 0 ) fWrappedProcess->SetProcessManager(mgr);
768 fWrappedProcess->SetProcessManager(mgr); << 807 else G4VProcess::SetProcessManager(mgr);
769 else <<
770 G4VProcess::SetProcessManager(mgr); <<
771 808
772 // -- initialize fSharedData pointer: 809 // -- initialize fSharedData pointer:
773 if (G4BiasingProcessSharedData::fSharedDataM << 810 if ( G4BiasingProcessSharedData::fSharedDataMap.Find(mgr) == G4BiasingProcessSharedData::fSharedDataMap.End() )
774 == G4BiasingProcessSharedData::fSharedDataM << 811 {
775 { << 812 fSharedData = new G4BiasingProcessSharedData( mgr );
776 fSharedData = new G4BiasingProcessSharedDa << 813 G4BiasingProcessSharedData::fSharedDataMap[mgr] = fSharedData;
777 G4BiasingProcessSharedData::fSharedDataMap << 814 }
778 } << 815 else fSharedData = G4BiasingProcessSharedData::fSharedDataMap[mgr] ;
779 else <<
780 { <<
781 fSharedData = G4BiasingProcessSharedData:: <<
782 } <<
783 // -- augment list of co-operating processes 816 // -- augment list of co-operating processes:
784 fSharedData->fBiasingProcessInterfaces.push_ << 817 fSharedData-> fBiasingProcessInterfaces.push_back( this );
785 fSharedData->fPublicBiasingProcessInterfaces << 818 fSharedData-> fPublicBiasingProcessInterfaces.push_back( this );
786 if ( fIsPhysicsBasedBiasing ) 819 if ( fIsPhysicsBasedBiasing )
787 { << 820 {
788 fSharedData->fPhysicsBiasingProcessInterfa << 821 fSharedData-> fPhysicsBiasingProcessInterfaces.push_back( this );
789 fSharedData-> fPublicPhysicsBiasingProcess << 822 fSharedData-> fPublicPhysicsBiasingProcessInterfaces.push_back( this );
790 } << 823 }
791 else 824 else
792 { << 825 {
793 fSharedData->fNonPhysicsBiasingProcessInte << 826 fSharedData-> fNonPhysicsBiasingProcessInterfaces.push_back( this );
794 fSharedData->fPublicNonPhysicsBiasingProce << 827 fSharedData-> fPublicNonPhysicsBiasingProcessInterfaces.push_back( this );
795 } << 828 }
796 // -- remember process manager: 829 // -- remember process manager:
797 fProcessManager = mgr; 830 fProcessManager = mgr;
798 } 831 }
799 832
>> 833
800 const G4ProcessManager* G4BiasingProcessInterf 834 const G4ProcessManager* G4BiasingProcessInterface::GetProcessManager()
801 { 835 {
802 if ( fWrappedProcess != nullptr ) << 836 if ( fWrappedProcess != 0 ) return fWrappedProcess->GetProcessManager();
803 return fWrappedProcess->GetProcessManager( << 837 else return G4VProcess::GetProcessManager();
804 else <<
805 return G4VProcess::GetProcessManager(); <<
806 } 838 }
807 839
808 void G4BiasingProcessInterface:: << 840
809 BuildWorkerPhysicsTable(const G4ParticleDefini << 841 void G4BiasingProcessInterface::BuildWorkerPhysicsTable(const G4ParticleDefinition& pd)
810 { 842 {
811 // -- Sequential mode : not called 843 // -- Sequential mode : not called
812 // -- MT mode : called after Prepare 844 // -- MT mode : called after PrepareWorkerPhysicsTable(..)
813 // -- PrepareWorkerPhysicsTable(...) has bee 845 // -- PrepareWorkerPhysicsTable(...) has been called first for all processes,
814 // -- so the first/last flags and G4BiasingP 846 // -- so the first/last flags and G4BiasingProcessInterface vector of processes have
815 // -- been properly setup, fIamFirstGPIL is 847 // -- been properly setup, fIamFirstGPIL is valid.
816 if ( fWrappedProcess != nullptr ) << 848 if ( fWrappedProcess != 0 )
817 { << 849 {
818 fWrappedProcess->BuildWorkerPhysicsTable(p << 850 fWrappedProcess->BuildWorkerPhysicsTable(pd);
819 } << 851 }
820 852
821 if ( fIamFirstGPIL ) 853 if ( fIamFirstGPIL )
822 { << 854 {
823 // -- Re-order vector of processes to matc << 855 // -- Re-order vector of processes to match that of the GPIL
824 // -- (made for fIamFirstGPIL, but importa << 856 // -- (made for fIamFirstGPIL, but important is to have it made once):
825 ReorderBiasingVectorAsGPIL(); << 857 ReorderBiasingVectorAsGPIL();
826 // -- Let operators to configure themselve << 858 // -- Let operators to configure themselves for the worker thread, if needed.
827 // -- Registration to physics model catalo << 859 // -- Registration to physics model catalog **IS NOT** to be made here, but in Configure().
828 // -- The fDoCommonConfigure is to ensure << 860 // -- The fDoCommonConfigure is to ensure that this Configure is made by only one process (othewise each first process makes the call):
829 if ( fDoCommonConfigure.Get() ) << 861 if ( fDoCommonConfigure.Get() )
830 { << 862 {
831 for ( std::size_t optr=0 ; optr<(G4VBias << 863 for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
832 { << 864 ( G4VBiasingOperator::GetBiasingOperators() )[optr]->ConfigureForWorker( );
833 (G4VBiasingOperator::GetBiasingOperato << 865 fDoCommonConfigure.Put(false);
834 } << 866 }
835 fDoCommonConfigure.Put(false); <<
836 } 867 }
837 } <<
838 } 868 }
839 869
840 void G4BiasingProcessInterface:: << 870
841 PrepareWorkerPhysicsTable(const G4ParticleDefi << 871 void G4BiasingProcessInterface::PrepareWorkerPhysicsTable(const G4ParticleDefinition& pd)
842 { 872 {
843 // -- Sequential mode : not called 873 // -- Sequential mode : not called
844 // -- MT mode : called first, before 874 // -- MT mode : called first, before BuildWorkerPhysicsTable(..)
845 // -- Let process finding its first/last pos 875 // -- Let process finding its first/last position in the process manager:
846 SetUpFirstLastFlags(); 876 SetUpFirstLastFlags();
847 877
848 if ( fWrappedProcess != nullptr ) << 878 if ( fWrappedProcess != 0 )
849 { << 879 {
850 fWrappedProcess->PrepareWorkerPhysicsTable << 880 fWrappedProcess->PrepareWorkerPhysicsTable(pd);
851 } << 881 }
852 } 882 }
853 883
>> 884
854 void G4BiasingProcessInterface::ResetNumberOfI 885 void G4BiasingProcessInterface::ResetNumberOfInteractionLengthLeft()
855 { 886 {
856 if ( fWrappedProcess != nullptr ) << 887 if ( fWrappedProcess != 0 ) fWrappedProcess->ResetNumberOfInteractionLengthLeft();
857 fWrappedProcess->ResetNumberOfInteractionL <<
858 } 888 }
859 889
860 G4bool G4BiasingProcessInterface:: << 890
861 GetIsFirstPostStepGPILInterface( G4bool physOn << 891 G4bool G4BiasingProcessInterface::GetIsFirstPostStepGPILInterface( G4bool physOnly ) const
862 { 892 {
863 G4int iPhys = ( physOnly ) ? 1 : 0; 893 G4int iPhys = ( physOnly ) ? 1 : 0;
864 return fFirstLastFlags[IdxFirstLast( 1, 1, i 894 return fFirstLastFlags[IdxFirstLast( 1, 1, iPhys)];
865 } 895 }
866 896
867 G4bool G4BiasingProcessInterface:: << 897
868 GetIsLastPostStepGPILInterface( G4bool physOnl << 898 G4bool G4BiasingProcessInterface::GetIsLastPostStepGPILInterface( G4bool physOnly ) const
869 { 899 {
870 G4int iPhys = ( physOnly ) ? 1 : 0; 900 G4int iPhys = ( physOnly ) ? 1 : 0;
871 return fFirstLastFlags[IdxFirstLast( 0, 1, i 901 return fFirstLastFlags[IdxFirstLast( 0, 1, iPhys)];
872 } 902 }
873 903
874 G4bool G4BiasingProcessInterface:: << 904
875 GetIsFirstPostStepDoItInterface( G4bool physOn << 905 G4bool G4BiasingProcessInterface::GetIsFirstPostStepDoItInterface( G4bool physOnly ) const
876 { 906 {
877 G4int iPhys = ( physOnly ) ? 1 : 0; 907 G4int iPhys = ( physOnly ) ? 1 : 0;
878 return fFirstLastFlags[IdxFirstLast( 1, 0, i 908 return fFirstLastFlags[IdxFirstLast( 1, 0, iPhys)];
879 } 909 }
880 910
881 G4bool G4BiasingProcessInterface:: << 911
882 GetIsLastPostStepDoItInterface( G4bool physOnl << 912 G4bool G4BiasingProcessInterface::GetIsLastPostStepDoItInterface( G4bool physOnly ) const
883 { 913 {
884 G4int iPhys = ( physOnly ) ? 1 : 0; 914 G4int iPhys = ( physOnly ) ? 1 : 0;
885 return fFirstLastFlags[IdxFirstLast( 0, 0, i 915 return fFirstLastFlags[IdxFirstLast( 0, 0, iPhys)];
886 } 916 }
887 917
888 G4bool G4BiasingProcessInterface:: << 918
889 IsFirstPostStepGPILInterface(G4bool physOnly) << 919 G4bool G4BiasingProcessInterface::IsFirstPostStepGPILInterface(G4bool physOnly) const
890 { 920 {
891 G4bool isFirst = true; 921 G4bool isFirst = true;
892 const G4ProcessVector* pv = fProcessManager- 922 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeGPIL);
893 G4int thisIdx(-1); 923 G4int thisIdx(-1);
894 for ( auto i = 0; i < (G4int)pv->size(); ++i << 924 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
895 { << 925 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
896 if ( (*pv)(i) == this ) { thisIdx = i; bre << 926 {
897 } << 927 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
898 if ( thisIdx < 0 ) return false; // -- to ig << 928 {
899 for ( std::size_t i=0; i<(fSharedData->fBias << 929 G4int thatIdx(-1);
900 { << 930 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
901 if ( (fSharedData->fBiasingProcessInterfac << 931 if ( thisIdx > thatIdx )
902 { << 932 {
903 G4int thatIdx(-1); << 933 isFirst = false;
904 for (auto j = 0; j < (G4int)pv->size(); << 934 break;
905 { << 935 }
906 if ( (*pv)(j) == (fSharedData->fBiasin << 936 }
907 { <<
908 thatIdx = j; break; <<
909 } <<
910 } <<
911 if ( thatIdx >= 0 ) // -- to ignore pure <<
912 { <<
913 if ( thisIdx > thatIdx ) <<
914 { <<
915 isFirst = false; <<
916 break; <<
917 } <<
918 } <<
919 } 937 }
920 } <<
921 return isFirst; 938 return isFirst;
922 } 939 }
923 940
924 G4bool G4BiasingProcessInterface:: << 941
925 IsLastPostStepGPILInterface(G4bool physOnly) c << 942 G4bool G4BiasingProcessInterface::IsLastPostStepGPILInterface(G4bool physOnly) const
926 { 943 {
927 G4bool isLast = true; 944 G4bool isLast = true;
928 const G4ProcessVector* pv = fProcessManager- 945 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeGPIL);
929 G4int thisIdx(-1); 946 G4int thisIdx(-1);
930 for (auto i = 0; i < (G4int)pv->size(); ++i << 947 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
931 { << 948 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
932 if ( (*pv)(i) == this ) { thisIdx = i; bre << 949 {
933 } << 950 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
934 if ( thisIdx < 0 ) return false; // -- to ig << 951 {
935 for (std::size_t i=0; i<(fSharedData->fBiasi << 952 G4int thatIdx(-1);
936 { << 953 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
937 if ( (fSharedData->fBiasingProcessInterfac << 954 if ( thisIdx < thatIdx )
938 { << 955 {
939 G4int thatIdx(-1); << 956 isLast = false;
940 for (auto j = 0; j < (G4int)pv->size(); << 957 break;
941 { << 958 }
942 if ( (*pv)(j) == (fSharedData->fBiasin << 959 }
943 { <<
944 thatIdx = j; break; <<
945 } <<
946 } <<
947 if ( thatIdx >= 0 ) // -- to ignore pure <<
948 { <<
949 if ( thisIdx < thatIdx ) <<
950 { <<
951 isLast = false; <<
952 break; <<
953 } <<
954 } <<
955 } 960 }
956 } <<
957 return isLast; 961 return isLast;
958 } 962 }
959 963
960 G4bool G4BiasingProcessInterface:: << 964
961 IsFirstPostStepDoItInterface(G4bool physOnly) << 965 G4bool G4BiasingProcessInterface::IsFirstPostStepDoItInterface(G4bool physOnly) const
962 { 966 {
963 G4bool isFirst = true; 967 G4bool isFirst = true;
964 const G4ProcessVector* pv = fProcessManager- 968 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeDoIt);
965 G4int thisIdx(-1); 969 G4int thisIdx(-1);
966 for (auto i = 0; i < (G4int)pv->size(); ++i << 970 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
967 { << 971 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
968 if ( (*pv)(i) == this ) { thisIdx = i; bre << 972 {
969 } << 973 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
970 if ( thisIdx < 0 ) return false; // -- to ig << 974 {
971 for (std::size_t i=0; i<(fSharedData->fBiasi << 975 G4int thatIdx(-1);
972 { << 976 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
973 if ( (fSharedData->fBiasingProcessInterfac << 977 if ( thisIdx > thatIdx )
974 { << 978 {
975 G4int thatIdx(-1); << 979 isFirst = false;
976 for (auto j = 0; j < (G4int)pv->size(); << 980 break;
977 { << 981 }
978 if ( (*pv)(j) == (fSharedData->fBiasin << 982 }
979 { <<
980 thatIdx = j; break; <<
981 } <<
982 } <<
983 if ( thatIdx >= 0 ) // -- to ignore pure <<
984 { <<
985 if ( thisIdx > thatIdx ) <<
986 { <<
987 isFirst = false; <<
988 break; <<
989 } <<
990 } <<
991 } 983 }
992 } <<
993 return isFirst; 984 return isFirst;
994 } 985 }
995 986
996 G4bool G4BiasingProcessInterface:: << 987
997 IsLastPostStepDoItInterface(G4bool physOnly) c << 988 G4bool G4BiasingProcessInterface::IsLastPostStepDoItInterface(G4bool physOnly) const
998 { 989 {
999 G4bool isLast = true; 990 G4bool isLast = true;
1000 const G4ProcessVector* pv = fProcessManager 991 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeDoIt);
1001 G4int thisIdx(-1); 992 G4int thisIdx(-1);
1002 for (auto i = 0; i < (G4int)pv->size(); ++i << 993 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
1003 { << 994 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
1004 if ( (*pv)(i) == this ) { thisIdx = i; br << 995 {
1005 } << 996 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
1006 if ( thisIdx < 0 ) return false; // -- to i << 997 {
1007 for (std::size_t i=0; i<(fSharedData->fBias << 998 G4int thatIdx(-1);
1008 { << 999 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
1009 if ( (fSharedData->fBiasingProcessInterfa << 1000 if ( thisIdx < thatIdx )
1010 { << 1001 {
1011 G4int thatIdx(-1); << 1002 isLast = false;
1012 for (auto j = 0; j < (G4int)pv->size(); << 1003 break;
1013 { << 1004 }
1014 if ( (*pv)(j) == (fSharedData->fBiasi << 1005 }
1015 { <<
1016 thatIdx = j; break; <<
1017 } <<
1018 } <<
1019 if ( thatIdx >= 0 ) // -- to ignore pur <<
1020 { <<
1021 if ( thisIdx < thatIdx ) <<
1022 { <<
1023 isLast = false; <<
1024 break; <<
1025 } <<
1026 } <<
1027 } 1006 }
1028 } <<
1029 return isLast; 1007 return isLast;
1030 } 1008 }
1031 1009
>> 1010
1032 void G4BiasingProcessInterface::SetUpFirstLas 1011 void G4BiasingProcessInterface::SetUpFirstLastFlags()
1033 { 1012 {
1034 for (G4int iPhys = 0; iPhys < 2; ++iPhys) << 1013 for ( G4int iPhys = 0; iPhys < 2; iPhys++ )
1035 { << 1014 {
1036 G4bool physOnly = ( iPhys == 1 ); << 1015 G4bool physOnly = ( iPhys == 1 );
1037 fFirstLastFlags[IdxFirstLast( 1, 1, iPhys << 1016 fFirstLastFlags[IdxFirstLast( 1, 1, iPhys)] = IsFirstPostStepGPILInterface(physOnly);
1038 fFirstLastFlags[IdxFirstLast( 0, 1, iPhys << 1017 fFirstLastFlags[IdxFirstLast( 0, 1, iPhys)] = IsLastPostStepGPILInterface(physOnly);
1039 fFirstLastFlags[IdxFirstLast( 1, 0, iPhys << 1018 fFirstLastFlags[IdxFirstLast( 1, 0, iPhys)] = IsFirstPostStepDoItInterface(physOnly);
1040 fFirstLastFlags[IdxFirstLast( 0, 0, iPhys << 1019 fFirstLastFlags[IdxFirstLast( 0, 0, iPhys)] = IsLastPostStepDoItInterface(physOnly);
1041 } << 1020 }
1042 1021
1043 // -- for itself, for optimization: 1022 // -- for itself, for optimization:
1044 fIamFirstGPIL = GetIsFirstPostStepGPILInte 1023 fIamFirstGPIL = GetIsFirstPostStepGPILInterface( false );
1045 } 1024 }
1046 1025
>> 1026
1047 void G4BiasingProcessInterface::ResetForUnbia 1027 void G4BiasingProcessInterface::ResetForUnbiasedTracking()
1048 { 1028 {
1049 fOccurenceBiasingOperation = nullptr; << 1029 fOccurenceBiasingOperation = 0;
1050 fFinalStateBiasingOperation = nullptr; << 1030 fFinalStateBiasingOperation = 0;
1051 fNonPhysicsBiasingOperation = nullptr; << 1031 fNonPhysicsBiasingOperation = 0;
1052 fBiasingInteractionLaw = nullptr; << 1032 fBiasingInteractionLaw = 0;
1053 } 1033 }
1054 1034
1055 void G4BiasingProcessInterface:: << 1035
1056 InvokeWrappedProcessPostStepGPIL( const G4Tra << 1036 void G4BiasingProcessInterface::InvokeWrappedProcessPostStepGPIL( const G4Track& track,
1057 G4double pr << 1037 G4double previousStepSize,
1058 G4ForceCond << 1038 G4ForceCondition* condition )
1059 { 1039 {
1060 G4double usedPreviousStepSize = previousSte 1040 G4double usedPreviousStepSize = previousStepSize;
1061 // -- if the physics process has been under << 1041 // -- if the physics process has been under occurence biasing in the previous step
1062 // -- we reset it, as we don't know if it w 1042 // -- we reset it, as we don't know if it will be biased again or not in this
1063 // -- step. The pity is that PostStepGPIL a 1043 // -- step. The pity is that PostStepGPIL and interaction length (cross-section)
1064 // -- calculations are done both in the Pos 1044 // -- calculations are done both in the PostStepGPIL of the process, while here we
1065 // -- are just interested in the calculatio 1045 // -- are just interested in the calculation of the cross-section. This is a pity
1066 // -- as this forces to re-generated a rand 1046 // -- as this forces to re-generated a random number for nothing.
1067 if ( fResetWrappedProcessInteractionLength 1047 if ( fResetWrappedProcessInteractionLength )
1068 { << 1048 {
1069 fResetWrappedProcessInteractionLength = f << 1049 fResetWrappedProcessInteractionLength = false;
1070 fWrappedProcess->ResetNumberOfInteraction << 1050 fWrappedProcess->ResetNumberOfInteractionLengthLeft();
1071 // -- We set "previous step size" as 0.0, << 1051 // -- We set "previous step size" as 0.0, to let the process believe this is first step:
1072 usedPreviousStepSize = 0.0; << 1052 usedPreviousStepSize = 0.0;
1073 } << 1053 }
1074 // -- GPIL response: 1054 // -- GPIL response:
1075 fWrappedProcessPostStepGPIL = fWrappedProce << 1055 fWrappedProcessPostStepGPIL = fWrappedProcess->PostStepGetPhysicalInteractionLength(track, usedPreviousStepSize, condition);
1076 fWrappedProcessForceCondition = *condition; << 1056 fWrappedProcessForceCondition = *condition;
1077 // -- and (inverse) cross-section: 1057 // -- and (inverse) cross-section:
1078 fWrappedProcessInteractionLength = fWrapped 1058 fWrappedProcessInteractionLength = fWrappedProcess->GetCurrentInteractionLength();
1079 } 1059 }
1080 1060
>> 1061
1081 void G4BiasingProcessInterface::ReorderBiasin 1062 void G4BiasingProcessInterface::ReorderBiasingVectorAsGPIL()
1082 { 1063 {
1083 // -- re-order vector of processes to match 1064 // -- re-order vector of processes to match that of the GPIL:
1084 std::vector < G4BiasingProcessInterface* > 1065 std::vector < G4BiasingProcessInterface* > tmpProcess ( fSharedData->fBiasingProcessInterfaces );
1085 ( fSharedData -> fBiasingProcessInterfaces 1066 ( fSharedData -> fBiasingProcessInterfaces ) . clear();
1086 ( fSharedData -> fPhysicsBiasingProcessInte 1067 ( fSharedData -> fPhysicsBiasingProcessInterfaces ) . clear();
1087 ( fSharedData -> fNonPhysicsBiasingProcessI 1068 ( fSharedData -> fNonPhysicsBiasingProcessInterfaces ) . clear();
1088 ( fSharedData -> fPublicBiasingProcessInter 1069 ( fSharedData -> fPublicBiasingProcessInterfaces ) . clear();
1089 ( fSharedData -> fPublicPhysicsBiasingProce 1070 ( fSharedData -> fPublicPhysicsBiasingProcessInterfaces ) . clear();
1090 ( fSharedData -> fPublicNonPhysicsBiasingPr 1071 ( fSharedData -> fPublicNonPhysicsBiasingProcessInterfaces ) . clear();
1091 1072
1092 const G4ProcessVector* pv = fProcessManager 1073 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeGPIL);
1093 for (auto i = 0; i < (G4int)pv->size(); ++i << 1074 for (G4int i = 0; i < pv->size(); i++ )
1094 { << 1075 {
1095 for (std::size_t j = 0; j < tmpProcess.si << 1076 for ( size_t j = 0; j < tmpProcess.size(); j++ )
1096 { << 1077 {
1097 if ( (*pv)(i) == tmpProcess[j] ) << 1078 if ( (*pv)(i) == tmpProcess[j] )
1098 { << 1079 {
1099 ( fSharedData->fBiasingProcessInterfa << 1080 ( fSharedData -> fBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1100 ( fSharedData->fPublicBiasingProcessI << 1081 ( fSharedData -> fPublicBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1101 if ( tmpProcess[j] -> fIsPhysicsBased << 1082 if ( tmpProcess[j] -> fIsPhysicsBasedBiasing )
1102 { << 1083 {
1103 ( fSharedData->fPhysicsBiasingProce << 1084 ( fSharedData -> fPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1104 ( fSharedData->fPublicPhysicsBiasin << 1085 ( fSharedData -> fPublicPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1105 } << 1086 }
1106 else << 1087 else
1107 { << 1088 {
1108 ( fSharedData -> fNonPhysicsBiasing << 1089 ( fSharedData -> fNonPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1109 ( fSharedData -> fPublicNonPhysicsB << 1090 ( fSharedData -> fPublicNonPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1110 } << 1091 }
1111 break; << 1092 break;
1112 } << 1093 }
>> 1094 }
1113 } 1095 }
1114 } <<
1115 } 1096 }
1116 1097