| 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 <<
350 { <<
351 *condition = NotForced; <<
352 return DBL_MAX; <<
353 } <<
354 } <<
355 391
356 // ----------------------------------------- 392 // --------------------------------------------------
357 // -- A biasing operator exists. Proceed wit 393 // -- A biasing operator exists. Proceed with
358 // -- treating non-physics and physics biasi 394 // -- treating non-physics and physics biasing cases:
359 //------------------------------------------ 395 //---------------------------------------------------
360 396
361 // -- non-physics-based biasing case: 397 // -- non-physics-based biasing case:
362 // ---------------------------------- 398 // ----------------------------------
363 if ( !fIsPhysicsBasedBiasing ) 399 if ( !fIsPhysicsBasedBiasing )
364 { << 400 {
365 fNonPhysicsBiasingOperation = (fSharedData << 401 fNonPhysicsBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedNonPhysicsBiasingOperation( &track, this );
366 if ( fNonPhysicsBiasingOperation == nullpt << 402 if ( fNonPhysicsBiasingOperation == 0 )
367 { << 403 {
368 *condition = NotForced; << 404 *condition = NotForced;
369 return DBL_MAX; << 405 return DBL_MAX;
>> 406 }
>> 407 return fNonPhysicsBiasingOperation->DistanceToApplyOperation(&track, previousStepSize, condition);
370 } 408 }
371 return fNonPhysicsBiasingOperation->Distan << 409
372 } <<
373 410
374 // -- Physics based biasing case: 411 // -- Physics based biasing case:
375 // ------------------------------ 412 // ------------------------------
376 // -- Ask for possible GPIL biasing operatio 413 // -- Ask for possible GPIL biasing operation:
377 fOccurenceBiasingOperation = (fSharedData->f 414 fOccurenceBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedOccurenceBiasingOperation( &track, this );
378 415
379 // -- no operation for occurrence biasing, a << 416
380 if ( fOccurenceBiasingOperation == nullptr ) << 417 // -- no operation for occurence biasing, analog GPIL returns the wrapped process GPIL and condition values
381 { << 418 if ( fOccurenceBiasingOperation == 0 )
382 *condition = fWrappedProcessForceCondition << 419 {
383 return fWrappedProcessPostStepGPIL; << 420 *condition = fWrappedProcessForceCondition;
384 } << 421 return fWrappedProcessPostStepGPIL;
>> 422 }
385 423
386 // -- A valid GPIL biasing operation has bee 424 // -- A valid GPIL biasing operation has been proposed:
387 // -- 0) remember wrapped process will need 425 // -- 0) remember wrapped process will need to be reset on biasing exit, if particle survives:
388 fResetWrappedProcessInteractionLength = true 426 fResetWrappedProcessInteractionLength = true;
389 // -- 1) update process interaction length f 427 // -- 1) update process interaction length for reference analog interaction law ( fWrappedProcessInteractionLength updated/collected above):
390 fPhysicalInteractionLaw->SetPhysicalCrossSec 428 fPhysicalInteractionLaw->SetPhysicalCrossSection( 1.0 / fWrappedProcessInteractionLength );
391 // -- 2) Collect biasing interaction law: 429 // -- 2) Collect biasing interaction law:
392 // -- The interaction law pointer is coll 430 // -- The interaction law pointer is collected as a const pointer to the interaction law object.
393 // -- This interaction law will be kept u 431 // -- This interaction law will be kept under control of the biasing operation, which is the only
394 // -- entity that will change the state o 432 // -- entity that will change the state of the biasing interaction law.
395 // -- The force condition for biasing is 433 // -- The force condition for biasing is asked at the same time, passing the analog one as default:
396 fBiasingForceCondition = fWrappedProcessForc << 434 fBiasingForceCondition = fWrappedProcessForceCondition;
397 fBiasingInteractionLaw = fOccurenceBiasingOp << 435 fBiasingInteractionLaw = fOccurenceBiasingOperation->ProvideOccurenceBiasingInteractionLaw( this, fBiasingForceCondition );
398 // -- 3) Ask operation to sample the biasing 436 // -- 3) Ask operation to sample the biasing interaction law:
399 fBiasingPostStepGPIL = fBiasingInteractionLa << 437 fBiasingPostStepGPIL = fBiasingInteractionLaw->GetSampledInteractionLength();
400 438
401 // -- finish 439 // -- finish
402 *condition = fBiasingForceCondition; 440 *condition = fBiasingForceCondition;
403 return fBiasingPostStepGPIL; << 441 return fBiasingPostStepGPIL;
>> 442
404 } 443 }
405 444
>> 445
>> 446
406 G4VParticleChange* G4BiasingProcessInterface:: 447 G4VParticleChange* G4BiasingProcessInterface::PostStepDoIt(const G4Track& track,
407 << 448 const G4Step& step)
408 { 449 {
409 // --------------------------------------- 450 // ---------------------------------------
410 // -- case outside of volume with biasing: 451 // -- case outside of volume with biasing:
411 // --------------------------------------- 452 // ---------------------------------------
412 if ( fSharedData->fCurrentBiasingOperator == << 453 if ( fSharedData->fCurrentBiasingOperator == 0 ) return fWrappedProcess->PostStepDoIt(track, step);
413 return fWrappedProcess->PostStepDoIt(track <<
414 454
415 // ---------------------------- 455 // ----------------------------
416 // -- non-physics biasing case: 456 // -- non-physics biasing case:
417 // ---------------------------- 457 // ----------------------------
418 if ( !fIsPhysicsBasedBiasing ) 458 if ( !fIsPhysicsBasedBiasing )
419 { << 459 {
420 G4VParticleChange* particleChange = fNonPh << 460 G4VParticleChange* particleChange = fNonPhysicsBiasingOperation->GenerateBiasingFinalState( &track, &step );
421 (fSharedData->fCurrentBiasingOperator)->Re << 461 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC_NonPhysics, fNonPhysicsBiasingOperation, particleChange );
422 return particleChange; << 462 return particleChange;
423 } << 463 }
424 464
425 // -- physics biasing case: 465 // -- physics biasing case:
426 // ------------------------ 466 // ------------------------
427 // -- It proceeds with the following logic: 467 // -- It proceeds with the following logic:
428 // -- 1) Obtain the final state 468 // -- 1) Obtain the final state
429 // -- This final state may be analog or b 469 // -- This final state may be analog or biased.
430 // -- The biased final state is obtained 470 // -- The biased final state is obtained through a biasing operator
431 // -- returned by the operator. 471 // -- returned by the operator.
432 // -- 2) The biased final state may be asked 472 // -- 2) The biased final state may be asked to be "force as it is"
433 // -- in what case the particle change is 473 // -- in what case the particle change is returned as is to the
434 // -- stepping. 474 // -- stepping.
435 // -- In all other cases (analog final st 475 // -- In all other cases (analog final state or biased final but
436 // -- not forced) the final state weight 476 // -- not forced) the final state weight may be modified by the
437 // -- occurrence biasing, if such an occu << 477 // -- occurence biasing, if such an occurence biasing is at play.
438 478
439 // -- Get final state, biased or analog: 479 // -- Get final state, biased or analog:
440 G4VParticleChange* finalStateParticleChange; << 480 G4VParticleChange* finalStateParticleChange;
441 G4BiasingAppliedCase BAC; 481 G4BiasingAppliedCase BAC;
442 fFinalStateBiasingOperation = (fSharedData-> 482 fFinalStateBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedFinalStateBiasingOperation( &track, this );
443 // -- Flag below is to force the biased gene << 483 // -- 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 << 484 // -- was or not a occurence biasing that would apply. Weight relevance under full responsibility of the biasing operation.
445 G4bool forceBiasedFinalState = false; 485 G4bool forceBiasedFinalState = false;
446 if ( fFinalStateBiasingOperation != nullptr << 486 if ( fFinalStateBiasingOperation != 0 )
447 { << 487 {
448 finalStateParticleChange = fFinalStateBias << 488 finalStateParticleChange = fFinalStateBiasingOperation->ApplyFinalStateBiasing( this, &track, &step, forceBiasedFinalState );
449 BAC = BAC_FinalState; << 489 BAC = BAC_FinalState;
450 } << 490 }
451 else 491 else
452 { << 492 {
453 finalStateParticleChange = fWrappedProcess << 493 finalStateParticleChange = fWrappedProcess->PostStepDoIt(track, step);
454 BAC = BAC_None ; << 494 BAC = BAC_None ;
455 } << 495 }
456 << 496
457 // -- if no occurrence biasing operation, we << 497 // -- if no occurence biasing operation, we're done:
458 if ( fOccurenceBiasingOperation == nullptr ) << 498 if ( fOccurenceBiasingOperation == 0 )
459 { << 499 {
460 (fSharedData->fCurrentBiasingOperator)->Re << 500 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC, fFinalStateBiasingOperation, finalStateParticleChange );
461 return finalStateParticleChange; << 501 return finalStateParticleChange;
462 } << 502 }
463 503
464 // -- if biased final state has been asked t 504 // -- if biased final state has been asked to be forced, we're done:
465 if ( forceBiasedFinalState ) 505 if ( forceBiasedFinalState )
466 { << 506 {
467 (fSharedData->fCurrentBiasingOperator)->Re << 507 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC, fFinalStateBiasingOperation, finalStateParticleChange );
468 return finalStateParticleChange; << 508 return finalStateParticleChange;
469 } << 509 }
>> 510
470 511
471 // -- If occurrence biasing, applies the occ << 512 // -- If occurence biasing, applies the occurence biasing weight correction on top of final state (biased or not):
472 G4double weightForInteraction = 1.0; 513 G4double weightForInteraction = 1.0;
473 if ( !fBiasingInteractionLaw->IsSingular() ) << 514 if ( !fBiasingInteractionLaw->IsSingular() ) weightForInteraction =
474 { << 515 fPhysicalInteractionLaw->ComputeEffectiveCrossSectionAt(step.GetStepLength()) /
475 weightForInteraction = fPhysicalInteractio << 516 fBiasingInteractionLaw ->ComputeEffectiveCrossSectionAt(step.GetStepLength());
476 / fBiasingInteraction <<
477 } <<
478 else 517 else
479 { << 518 {
480 // -- at this point effective XS can only << 519 // -- at this point effective XS can only be infinite, if not, there is a logic problem
481 if ( !fBiasingInteractionLaw->IsEffectiveC << 520 if ( !fBiasingInteractionLaw->IsEffectiveCrossSectionInfinite() )
>> 521 {
>> 522 G4ExceptionDescription ed;
>> 523 ed << "Internal inconsistency in cross-section handling. Please report !" << G4endl;
>> 524 G4Exception(" G4BiasingProcessInterface::PostStepDoIt(...)",
>> 525 "BIAS.GEN.02",
>> 526 JustWarning,
>> 527 ed);
>> 528 // -- if XS is infinite, weight is zero (and will stay zero), but we'll do differently.
>> 529 // -- Should foresee in addition something to remember that in case of singular
>> 530 // -- distribution, weight can only be partly calculated
>> 531 }
>> 532 }
>> 533
>> 534 if ( weightForInteraction <= 0. )
482 { 535 {
483 G4ExceptionDescription ed; 536 G4ExceptionDescription ed;
484 ed << "Internal inconsistency in cross-s << 537 ed << " Negative interaction weight : w_I = "
>> 538 << weightForInteraction <<
>> 539 " XS_I(phys) = " << fBiasingInteractionLaw ->ComputeEffectiveCrossSectionAt(step.GetStepLength()) <<
>> 540 " XS_I(bias) = " << fPhysicalInteractionLaw->ComputeEffectiveCrossSectionAt(step.GetStepLength()) <<
>> 541 " step length = " << step.GetStepLength() <<
>> 542 " Interaction law = `" << fBiasingInteractionLaw << "'" <<
>> 543 G4endl;
485 G4Exception(" G4BiasingProcessInterface: 544 G4Exception(" G4BiasingProcessInterface::PostStepDoIt(...)",
486 "BIAS.GEN.02", JustWarning, << 545 "BIAS.GEN.03",
487 // -- if XS is infinite, weight is zero << 546 JustWarning,
488 // -- Should foresee in addition somethi << 547 ed);
489 // -- distribution, weight can only be p << 548
490 } 549 }
491 } <<
492 550
493 if ( weightForInteraction <= 0. ) << 551 (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC,
494 { << 552 fOccurenceBiasingOperation, weightForInteraction,
495 G4ExceptionDescription ed; << 553 fFinalStateBiasingOperation, finalStateParticleChange );
496 ed << " Negative interaction weight : w_I << 554
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 555
514 fOccurenceBiasingParticleChange->SetOccurenc 556 fOccurenceBiasingParticleChange->SetOccurenceWeightForInteraction( weightForInteraction );
515 fOccurenceBiasingParticleChange->SetSecondar 557 fOccurenceBiasingParticleChange->SetSecondaryWeightByProcess( true );
516 fOccurenceBiasingParticleChange->SetWrappedP 558 fOccurenceBiasingParticleChange->SetWrappedParticleChange( finalStateParticleChange );
517 fOccurenceBiasingParticleChange->ProposeTrac 559 fOccurenceBiasingParticleChange->ProposeTrackStatus( finalStateParticleChange->GetTrackStatus() );
518 fOccurenceBiasingParticleChange->StealSecond 560 fOccurenceBiasingParticleChange->StealSecondaries(); // -- this also makes weightForInteraction applied to secondaries stolen
519 561
520 // -- finish: 562 // -- finish:
521 return fOccurenceBiasingParticleChange; 563 return fOccurenceBiasingParticleChange;
>> 564
522 } 565 }
523 566
>> 567
524 // -- AlongStep methods: 568 // -- AlongStep methods:
525 G4double G4BiasingProcessInterface:: << 569 G4double G4BiasingProcessInterface::AlongStepGetPhysicalInteractionLength(const G4Track& track,
526 AlongStepGetPhysicalInteractionLength(const G4 << 570 G4double previousStepSize,
527 G4 << 571 G4double currentMinimumStep,
528 G4 << 572 G4double& proposedSafety,
529 G4 << 573 G4GPILSelection* selection)
530 G4 <<
531 { 574 {
>> 575
532 // -- for helper methods: 576 // -- for helper methods:
533 fCurrentMinimumStep = currentMinimumStep; 577 fCurrentMinimumStep = currentMinimumStep;
534 fProposedSafety = proposedSafety; << 578 fProposedSafety = proposedSafety;
>> 579
535 580
536 // -- initialization default case: 581 // -- initialization default case:
537 fWrappedProcessAlongStepGPIL = DBL_MAX; << 582 fWrappedProcessAlongStepGPIL = DBL_MAX;
538 *selection = NotCandidateF 583 *selection = NotCandidateForSelection;
539 // --------------------------------------- 584 // ---------------------------------------
540 // -- case outside of volume with biasing: 585 // -- case outside of volume with biasing:
541 // --------------------------------------- 586 // ---------------------------------------
542 if ( fSharedData->fCurrentBiasingOperator == << 587 if ( fSharedData->fCurrentBiasingOperator == 0 )
543 { << 588 {
544 if ( fWrappedProcessIsAlong ) << 589 if ( fWrappedProcessIsAlong ) fWrappedProcessAlongStepGPIL =
545 fWrappedProcessAlongStepGPIL = fWrappedP << 590 fWrappedProcess->AlongStepGetPhysicalInteractionLength(track,
546 ->AlongStepGetPhysicalInteractionLengt << 591 previousStepSize,
547 << 592 currentMinimumStep,
548 << 593 proposedSafety,
549 return fWrappedProcessAlongStepGPIL; << 594 selection);
550 } << 595 return fWrappedProcessAlongStepGPIL;
>> 596 }
551 597
552 // ----------------------------------------- 598 // --------------------------------------------------------------------
553 // -- non-physics based biasing: no along op 599 // -- non-physics based biasing: no along operation expected (for now):
554 // ----------------------------------------- 600 // --------------------------------------------------------------------
555 if ( !fIsPhysicsBasedBiasing ) return fWrapp 601 if ( !fIsPhysicsBasedBiasing ) return fWrappedProcessAlongStepGPIL;
556 602
557 // ---------------------- 603 // ----------------------
558 // -- physics-based case: 604 // -- physics-based case:
559 // ---------------------- 605 // ----------------------
560 if ( fOccurenceBiasingOperation == nullptr ) << 606 if ( fOccurenceBiasingOperation == 0 )
561 { << 607 {
562 if ( fWrappedProcessIsAlong ) << 608 if ( fWrappedProcessIsAlong ) fWrappedProcessAlongStepGPIL =
563 fWrappedProcessAlongStepGPIL = fWrappedP << 609 fWrappedProcess->AlongStepGetPhysicalInteractionLength(track,
564 ->AlongStepGetPhysicalInteractionLengt << 610 previousStepSize,
565 << 611 currentMinimumStep,
566 << 612 proposedSafety,
567 return fWrappedProcessAlongStepGPIL; << 613 selection);
568 } << 614 return fWrappedProcessAlongStepGPIL;
569 << 615 }
570 // ----------------------------------------- << 616
571 // -- From here we have an valid occurrence << 617
572 // ----------------------------------------- << 618 // ----------------------------------------------------------
>> 619 // -- From here we have an valid occurence biasing operation:
>> 620 // ----------------------------------------------------------
573 // -- Give operation opportunity to shorten 621 // -- Give operation opportunity to shorten step proposed by physics process:
574 fBiasingAlongStepGPIL = fOccurenceBiasingOpe << 622 fBiasingAlongStepGPIL = fOccurenceBiasingOperation->ProposeAlongStepLimit( this );
575 G4double minimumStep = fBiasingAlongStepGPI << 623 G4double minimumStep = fBiasingAlongStepGPIL < currentMinimumStep ? fBiasingAlongStepGPIL : currentMinimumStep ;
576 ? fBiasingAlongStepGPI <<
577 // -- wrapped process is called with minimum 624 // -- wrapped process is called with minimum step ( <= currentMinimumStep passed ) : an along process can not
578 // -- have its operation stretched over what 625 // -- have its operation stretched over what it expects:
579 if ( fWrappedProcessIsAlong ) 626 if ( fWrappedProcessIsAlong )
580 { << 627 {
581 fWrappedProcessAlongStepGPIL = fWrappedPro << 628 fWrappedProcessAlongStepGPIL = fWrappedProcess->AlongStepGetPhysicalInteractionLength(track,
582 ->AlongStepGetPhysicalInteractionLength( << 629 previousStepSize,
583 << 630 minimumStep,
584 << 631 proposedSafety,
585 fWrappedProcessGPILSelection = *selection; << 632 selection);
586 fBiasingGPILSelection = fOccurenceBiasingO << 633 fWrappedProcessGPILSelection = *selection;
587 ->ProposeGPILSelection( fWrappedProcessG << 634 fBiasingGPILSelection = fOccurenceBiasingOperation->ProposeGPILSelection( fWrappedProcessGPILSelection );
588 } << 635 }
589 else 636 else
590 { << 637 {
591 fBiasingGPILSelection = fOccurenceBiasingO << 638 fBiasingGPILSelection = fOccurenceBiasingOperation->ProposeGPILSelection( NotCandidateForSelection );
592 ->ProposeGPILSelection( NotCandidateForS << 639 fWrappedProcessAlongStepGPIL = fBiasingAlongStepGPIL;
593 fWrappedProcessAlongStepGPIL = fBiasingAlo << 640 }
594 } <<
595 641
596 *selection = fBiasingGPILSelection; 642 *selection = fBiasingGPILSelection;
597 643
598 return fWrappedProcessAlongStepGPIL; 644 return fWrappedProcessAlongStepGPIL;
>> 645
599 } 646 }
600 647
601 G4VParticleChange* << 648 G4VParticleChange* G4BiasingProcessInterface::AlongStepDoIt(const G4Track& track,
602 G4BiasingProcessInterface::AlongStepDoIt(const << 649 const G4Step& step)
603 const <<
604 { 650 {
>> 651
605 // --------------------------------------- 652 // ---------------------------------------
606 // -- case outside of volume with biasing: 653 // -- case outside of volume with biasing:
607 // --------------------------------------- 654 // ---------------------------------------
608 if ( fSharedData->fCurrentBiasingOperator == << 655 if ( fSharedData->fCurrentBiasingOperator == 0 )
609 { <<
610 if ( fWrappedProcessIsAlong ) <<
611 { 656 {
612 return fWrappedProcess->AlongStepDoIt(tr << 657 if ( fWrappedProcessIsAlong ) return fWrappedProcess->AlongStepDoIt(track, step);
>> 658 else
>> 659 {
>> 660 fDummyParticleChange->Initialize( track );
>> 661 return fDummyParticleChange;
>> 662 }
613 } 663 }
614 else <<
615 { <<
616 fDummyParticleChange->Initialize( track <<
617 return fDummyParticleChange; <<
618 } <<
619 } <<
620 664
621 // ----------------------------------- 665 // -----------------------------------
622 // -- case inside volume with biasing: 666 // -- case inside volume with biasing:
623 // ----------------------------------- 667 // -----------------------------------
624 if ( fWrappedProcessIsAlong ) << 668 if ( fWrappedProcessIsAlong ) fOccurenceBiasingParticleChange->SetWrappedParticleChange( fWrappedProcess->AlongStepDoIt(track, step) );
625 { <<
626 fOccurenceBiasingParticleChange <<
627 ->SetWrappedParticleChange(fWrappedProce <<
628 } <<
629 else 669 else
630 { << 670 {
631 fOccurenceBiasingParticleChange->SetWrappe << 671 fOccurenceBiasingParticleChange->SetWrappedParticleChange ( 0 );
632 fOccurenceBiasingParticleChange->ProposeTr << 672 fOccurenceBiasingParticleChange->ProposeTrackStatus( track.GetTrackStatus() );
633 } << 673 }
634 G4double weightForNonInteraction (1.0); 674 G4double weightForNonInteraction (1.0);
635 if ( fBiasingInteractionLaw != nullptr ) << 675 if ( fBiasingInteractionLaw != 0 )
636 { << 676 {
637 weightForNonInteraction = << 677 weightForNonInteraction =
638 fPhysicalInteractionLaw->ComputeNonInter << 678 fPhysicalInteractionLaw->ComputeNonInteractionProbabilityAt(step.GetStepLength()) /
639 fBiasingInteractionLaw ->ComputeNonInter << 679 fBiasingInteractionLaw ->ComputeNonInteractionProbabilityAt(step.GetStepLength());
640 680
641 fOccurenceBiasingOperation->AlongMoveBy( t << 681 fOccurenceBiasingOperation->AlongMoveBy( this, &step, weightForNonInteraction );
642 682
643 if ( weightForNonInteraction <= 0. ) << 683 if ( weightForNonInteraction <= 0. )
644 { << 684 {
645 G4ExceptionDescription ed; << 685 G4ExceptionDescription ed;
646 ed << " Negative non interaction weight << 686 ed << " Negative non interaction weight : w_NI = " << weightForNonInteraction <<
647 " p_NI(phys) = " << fPhysicalInte << 687 " p_NI(phys) = " << fPhysicalInteractionLaw->ComputeNonInteractionProbabilityAt(step.GetStepLength()) <<
648 " p_NI(bias) = " << fBiasingInter << 688 " p_NI(bias) = " << fBiasingInteractionLaw ->ComputeNonInteractionProbabilityAt(step.GetStepLength()) <<
649 " step length = " << step.GetSte << 689 " step length = " << step.GetStepLength() <<
650 " biasing interaction law = `" << << 690 " biasing interaction law = `" << fBiasingInteractionLaw->GetName() << "'" << G4endl;
651 G4Exception(" G4BiasingProcessInterface: << 691 G4Exception(" G4BiasingProcessInterface::AlongStepDoIt(...)",
652 "BIAS.GEN.04", JustWarning, << 692 "BIAS.GEN.04",
>> 693 JustWarning,
>> 694 ed);
>> 695 }
>> 696
653 } 697 }
654 } <<
655 698
656 fOccurenceBiasingParticleChange << 699 fOccurenceBiasingParticleChange->SetOccurenceWeightForNonInteraction( weightForNonInteraction );
657 ->SetOccurenceWeightForNonInteraction( wei <<
658 700
659 return fOccurenceBiasingParticleChange; 701 return fOccurenceBiasingParticleChange;
>> 702
660 } 703 }
661 704
662 // -- AtRest methods 705 // -- AtRest methods
663 G4double G4BiasingProcessInterface:: << 706 G4double G4BiasingProcessInterface::AtRestGetPhysicalInteractionLength(const G4Track& track,
664 AtRestGetPhysicalInteractionLength(const G4Tra << 707 G4ForceCondition* condition)
665 G4ForceCond <<
666 { 708 {
667 return fWrappedProcess->AtRestGetPhysicalInt << 709 return fWrappedProcess->AtRestGetPhysicalInteractionLength(track, condition);
668 } 710 }
669 <<
670 G4VParticleChange* G4BiasingProcessInterface:: 711 G4VParticleChange* G4BiasingProcessInterface::AtRestDoIt(const G4Track& track,
671 << 712 const G4Step& step)
672 { 713 {
673 return fWrappedProcess->AtRestDoIt(track, st << 714 return fWrappedProcess->AtRestDoIt(track, step);
674 } 715 }
675 716
676 G4bool G4BiasingProcessInterface::IsApplicable << 717
>> 718 G4bool G4BiasingProcessInterface::IsApplicable(const G4ParticleDefinition& pd)
677 { 719 {
678 if ( fWrappedProcess != nullptr ) << 720 if ( fWrappedProcess != 0 ) return fWrappedProcess->IsApplicable(pd);
679 return fWrappedProcess->IsApplicable(pd); << 721 else return true;
680 else <<
681 return true; <<
682 } 722 }
683 723
684 void G4BiasingProcessInterface::SetMasterProc << 724
>> 725 void G4BiasingProcessInterface::SetMasterProcess(G4VProcess* masterP)
685 { 726 {
686 // -- Master for this process: 727 // -- Master for this process:
687 G4VProcess::SetMasterProcess(masterP); 728 G4VProcess::SetMasterProcess(masterP);
688 // -- Master for wrapped process: 729 // -- Master for wrapped process:
689 if ( fWrappedProcess != nullptr ) << 730 if ( fWrappedProcess != 0 )
690 { << 731 {
691 const G4BiasingProcessInterface* thisWrapp << 732 const G4BiasingProcessInterface* thisWrapperMaster = (const G4BiasingProcessInterface *)GetMasterProcess();
692 = (const G4BiasingProcessInterface *)Get << 733 // -- paranoia check: (?)
693 // -- paranoia check: (?) << 734 G4VProcess* wrappedMaster = 0;
694 G4VProcess* wrappedMaster = nullptr; << 735 wrappedMaster = thisWrapperMaster->GetWrappedProcess();
695 wrappedMaster = thisWrapperMaster->GetWrap << 736 fWrappedProcess->SetMasterProcess( wrappedMaster );
696 fWrappedProcess->SetMasterProcess( wrapped << 737 }
697 } <<
698 } 738 }
699 739
700 void G4BiasingProcessInterface:: << 740
701 BuildPhysicsTable(const G4ParticleDefinition& << 741 void G4BiasingProcessInterface::BuildPhysicsTable(const G4ParticleDefinition& pd)
702 { 742 {
703 // -- Sequential mode : called second (after 743 // -- Sequential mode : called second (after PreparePhysicsTable(..))
704 // -- MT mode : called second (after 744 // -- MT mode : called second (after PreparePhysicsTable(..)) by master thread.
705 // -- Corresponding proces 745 // -- Corresponding process instance not used then by tracking.
706 // -- PreparePhysicsTable(...) has been call 746 // -- PreparePhysicsTable(...) has been called first for all processes,
707 // -- so the first/last flags and G4BiasingP 747 // -- so the first/last flags and G4BiasingProcessInterface vector of processes have
708 // -- been properly setup, fIamFirstGPIL is 748 // -- been properly setup, fIamFirstGPIL is valid.
709 if ( fWrappedProcess != nullptr ) << 749 if ( fWrappedProcess != 0 )
710 { << 750 {
711 fWrappedProcess->BuildPhysicsTable(pd); << 751 fWrappedProcess->BuildPhysicsTable(pd);
712 } << 752 }
713 753
714 if ( fIamFirstGPIL ) 754 if ( fIamFirstGPIL )
715 { << 755 {
716 // -- Re-order vector of processes to matc << 756 // -- Re-order vector of processes to match that of the GPIL
717 // -- (made for fIamFirstGPIL, but importa << 757 // -- (made for fIamFirstGPIL, but important is to have it made once):
718 ReorderBiasingVectorAsGPIL(); << 758 ReorderBiasingVectorAsGPIL();
719 // -- Let operators to configure themselve << 759 // -- Let operators to configure themselves for the master thread or for sequential mode.
720 // -- Intended here is in particular the r << 760 // -- Intended here is in particular the registration to physics model catalog.
721 // -- The fDoCommonConfigure is to ensure << 761 // -- 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() ) << 762 if ( fDoCommonConfigure.Get() )
723 { << 763 {
724 for ( std::size_t optr=0; optr<(G4VBiasi << 764 for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
725 { << 765 ( G4VBiasingOperator::GetBiasingOperators() )[optr]->Configure( );
726 (G4VBiasingOperator::GetBiasingOperato << 766 fDoCommonConfigure.Put(false);
727 } << 767 }
728 fDoCommonConfigure.Put(false); << 768
729 } 769 }
730 } <<
731 } 770 }
732 771
733 void G4BiasingProcessInterface:: << 772
734 PreparePhysicsTable(const G4ParticleDefinition << 773 void G4BiasingProcessInterface::PreparePhysicsTable(const G4ParticleDefinition& pd)
735 { 774 {
736 // -- Sequential mode : called first (before 775 // -- Sequential mode : called first (before BuildPhysicsTable(..))
737 // -- MT mode : called first (before 776 // -- MT mode : called first (before BuildPhysicsTable(..)) by master thread.
738 // -- Corresponding proces 777 // -- Corresponding process instance not used then by tracking.
739 // -- Let process finding its first/last pos 778 // -- Let process finding its first/last position in the process manager:
740 SetUpFirstLastFlags(); 779 SetUpFirstLastFlags();
741 if ( fWrappedProcess != nullptr ) << 780 if ( fWrappedProcess != 0 )
742 { << 781 {
743 fWrappedProcess->PreparePhysicsTable(pd); << 782 fWrappedProcess->PreparePhysicsTable(pd);
744 } << 783 }
745 } 784 }
746 785
747 G4bool G4BiasingProcessInterface:: << 786
748 StorePhysicsTable(const G4ParticleDefinition* << 787 G4bool G4BiasingProcessInterface::StorePhysicsTable(const G4ParticleDefinition* pd, const G4String& s, G4bool f)
749 { 788 {
750 if ( fWrappedProcess != nullptr ) << 789 if ( fWrappedProcess != 0 ) return fWrappedProcess->StorePhysicsTable(pd, s, f);
751 return fWrappedProcess->StorePhysicsTable( << 790 else return false;
752 else <<
753 return false; <<
754 } 791 }
755 792
756 G4bool G4BiasingProcessInterface:: << 793
757 RetrievePhysicsTable(const G4ParticleDefinitio << 794 G4bool G4BiasingProcessInterface::RetrievePhysicsTable(const G4ParticleDefinition* pd, const G4String& s, G4bool f)
758 { 795 {
759 if ( fWrappedProcess != nullptr ) << 796 if ( fWrappedProcess != 0 ) return fWrappedProcess->RetrievePhysicsTable(pd, s, f);
760 return fWrappedProcess->RetrievePhysicsTab << 797 else return false;
761 else <<
762 return false; <<
763 } 798 }
764 799
>> 800
765 void G4BiasingProcessInterface::SetProcessMana 801 void G4BiasingProcessInterface::SetProcessManager(const G4ProcessManager* mgr)
766 { 802 {
767 if ( fWrappedProcess != nullptr ) << 803 if ( fWrappedProcess != 0 ) fWrappedProcess->SetProcessManager(mgr);
768 fWrappedProcess->SetProcessManager(mgr); << 804 else G4VProcess::SetProcessManager(mgr);
769 else <<
770 G4VProcess::SetProcessManager(mgr); <<
771 805
772 // -- initialize fSharedData pointer: 806 // -- initialize fSharedData pointer:
773 if (G4BiasingProcessSharedData::fSharedDataM << 807 if ( G4BiasingProcessSharedData::fSharedDataMap.Find(mgr) == G4BiasingProcessSharedData::fSharedDataMap.End() )
774 == G4BiasingProcessSharedData::fSharedDataM << 808 {
775 { << 809 fSharedData = new G4BiasingProcessSharedData( mgr );
776 fSharedData = new G4BiasingProcessSharedDa << 810 G4BiasingProcessSharedData::fSharedDataMap[mgr] = fSharedData;
777 G4BiasingProcessSharedData::fSharedDataMap << 811 }
778 } << 812 else fSharedData = G4BiasingProcessSharedData::fSharedDataMap[mgr] ;
779 else <<
780 { <<
781 fSharedData = G4BiasingProcessSharedData:: <<
782 } <<
783 // -- augment list of co-operating processes 813 // -- augment list of co-operating processes:
784 fSharedData->fBiasingProcessInterfaces.push_ << 814 fSharedData-> fBiasingProcessInterfaces.push_back( this );
785 fSharedData->fPublicBiasingProcessInterfaces << 815 fSharedData-> fPublicBiasingProcessInterfaces.push_back( this );
786 if ( fIsPhysicsBasedBiasing ) 816 if ( fIsPhysicsBasedBiasing )
787 { << 817 {
788 fSharedData->fPhysicsBiasingProcessInterfa << 818 fSharedData-> fPhysicsBiasingProcessInterfaces.push_back( this );
789 fSharedData-> fPublicPhysicsBiasingProcess << 819 fSharedData-> fPublicPhysicsBiasingProcessInterfaces.push_back( this );
790 } << 820 }
791 else 821 else
792 { << 822 {
793 fSharedData->fNonPhysicsBiasingProcessInte << 823 fSharedData-> fNonPhysicsBiasingProcessInterfaces.push_back( this );
794 fSharedData->fPublicNonPhysicsBiasingProce << 824 fSharedData-> fPublicNonPhysicsBiasingProcessInterfaces.push_back( this );
795 } << 825 }
796 // -- remember process manager: 826 // -- remember process manager:
797 fProcessManager = mgr; 827 fProcessManager = mgr;
798 } 828 }
799 829
>> 830
800 const G4ProcessManager* G4BiasingProcessInterf 831 const G4ProcessManager* G4BiasingProcessInterface::GetProcessManager()
801 { 832 {
802 if ( fWrappedProcess != nullptr ) << 833 if ( fWrappedProcess != 0 ) return fWrappedProcess->GetProcessManager();
803 return fWrappedProcess->GetProcessManager( << 834 else return G4VProcess::GetProcessManager();
804 else <<
805 return G4VProcess::GetProcessManager(); <<
806 } 835 }
807 836
808 void G4BiasingProcessInterface:: << 837
809 BuildWorkerPhysicsTable(const G4ParticleDefini << 838 void G4BiasingProcessInterface::BuildWorkerPhysicsTable(const G4ParticleDefinition& pd)
810 { 839 {
811 // -- Sequential mode : not called 840 // -- Sequential mode : not called
812 // -- MT mode : called after Prepare 841 // -- MT mode : called after PrepareWorkerPhysicsTable(..)
813 // -- PrepareWorkerPhysicsTable(...) has bee 842 // -- PrepareWorkerPhysicsTable(...) has been called first for all processes,
814 // -- so the first/last flags and G4BiasingP 843 // -- so the first/last flags and G4BiasingProcessInterface vector of processes have
815 // -- been properly setup, fIamFirstGPIL is 844 // -- been properly setup, fIamFirstGPIL is valid.
816 if ( fWrappedProcess != nullptr ) << 845 if ( fWrappedProcess != 0 )
817 { << 846 {
818 fWrappedProcess->BuildWorkerPhysicsTable(p << 847 fWrappedProcess->BuildWorkerPhysicsTable(pd);
819 } << 848 }
820 849
821 if ( fIamFirstGPIL ) 850 if ( fIamFirstGPIL )
822 { << 851 {
823 // -- Re-order vector of processes to matc << 852 // -- Re-order vector of processes to match that of the GPIL
824 // -- (made for fIamFirstGPIL, but importa << 853 // -- (made for fIamFirstGPIL, but important is to have it made once):
825 ReorderBiasingVectorAsGPIL(); << 854 ReorderBiasingVectorAsGPIL();
826 // -- Let operators to configure themselve << 855 // -- Let operators to configure themselves for the worker thread, if needed.
827 // -- Registration to physics model catalo << 856 // -- Registration to physics model catalog **IS NOT** to be made here, but in Configure().
828 // -- The fDoCommonConfigure is to ensure << 857 // -- 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() ) << 858 if ( fDoCommonConfigure.Get() )
830 { << 859 {
831 for ( std::size_t optr=0 ; optr<(G4VBias << 860 for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
832 { << 861 ( G4VBiasingOperator::GetBiasingOperators() )[optr]->ConfigureForWorker( );
833 (G4VBiasingOperator::GetBiasingOperato << 862 fDoCommonConfigure.Put(false);
834 } << 863 }
835 fDoCommonConfigure.Put(false); <<
836 } 864 }
837 } <<
838 } 865 }
839 866
840 void G4BiasingProcessInterface:: << 867
841 PrepareWorkerPhysicsTable(const G4ParticleDefi << 868 void G4BiasingProcessInterface::PrepareWorkerPhysicsTable(const G4ParticleDefinition& pd)
842 { 869 {
843 // -- Sequential mode : not called 870 // -- Sequential mode : not called
844 // -- MT mode : called first, before 871 // -- MT mode : called first, before BuildWorkerPhysicsTable(..)
845 // -- Let process finding its first/last pos 872 // -- Let process finding its first/last position in the process manager:
846 SetUpFirstLastFlags(); 873 SetUpFirstLastFlags();
847 874
848 if ( fWrappedProcess != nullptr ) << 875 if ( fWrappedProcess != 0 )
849 { << 876 {
850 fWrappedProcess->PrepareWorkerPhysicsTable << 877 fWrappedProcess->PrepareWorkerPhysicsTable(pd);
851 } << 878 }
852 } 879 }
853 880
>> 881
854 void G4BiasingProcessInterface::ResetNumberOfI 882 void G4BiasingProcessInterface::ResetNumberOfInteractionLengthLeft()
855 { 883 {
856 if ( fWrappedProcess != nullptr ) << 884 if ( fWrappedProcess != 0 ) fWrappedProcess->ResetNumberOfInteractionLengthLeft();
857 fWrappedProcess->ResetNumberOfInteractionL <<
858 } 885 }
859 886
860 G4bool G4BiasingProcessInterface:: << 887
861 GetIsFirstPostStepGPILInterface( G4bool physOn << 888 G4bool G4BiasingProcessInterface::GetIsFirstPostStepGPILInterface( G4bool physOnly ) const
862 { 889 {
863 G4int iPhys = ( physOnly ) ? 1 : 0; 890 G4int iPhys = ( physOnly ) ? 1 : 0;
864 return fFirstLastFlags[IdxFirstLast( 1, 1, i 891 return fFirstLastFlags[IdxFirstLast( 1, 1, iPhys)];
865 } 892 }
866 893
867 G4bool G4BiasingProcessInterface:: << 894
868 GetIsLastPostStepGPILInterface( G4bool physOnl << 895 G4bool G4BiasingProcessInterface::GetIsLastPostStepGPILInterface( G4bool physOnly ) const
869 { 896 {
870 G4int iPhys = ( physOnly ) ? 1 : 0; 897 G4int iPhys = ( physOnly ) ? 1 : 0;
871 return fFirstLastFlags[IdxFirstLast( 0, 1, i 898 return fFirstLastFlags[IdxFirstLast( 0, 1, iPhys)];
872 } 899 }
873 900
874 G4bool G4BiasingProcessInterface:: << 901
875 GetIsFirstPostStepDoItInterface( G4bool physOn << 902 G4bool G4BiasingProcessInterface::GetIsFirstPostStepDoItInterface( G4bool physOnly ) const
876 { 903 {
877 G4int iPhys = ( physOnly ) ? 1 : 0; 904 G4int iPhys = ( physOnly ) ? 1 : 0;
878 return fFirstLastFlags[IdxFirstLast( 1, 0, i 905 return fFirstLastFlags[IdxFirstLast( 1, 0, iPhys)];
879 } 906 }
880 907
881 G4bool G4BiasingProcessInterface:: << 908
882 GetIsLastPostStepDoItInterface( G4bool physOnl << 909 G4bool G4BiasingProcessInterface::GetIsLastPostStepDoItInterface( G4bool physOnly ) const
883 { 910 {
884 G4int iPhys = ( physOnly ) ? 1 : 0; 911 G4int iPhys = ( physOnly ) ? 1 : 0;
885 return fFirstLastFlags[IdxFirstLast( 0, 0, i 912 return fFirstLastFlags[IdxFirstLast( 0, 0, iPhys)];
886 } 913 }
887 914
888 G4bool G4BiasingProcessInterface:: << 915
889 IsFirstPostStepGPILInterface(G4bool physOnly) << 916 G4bool G4BiasingProcessInterface::IsFirstPostStepGPILInterface(G4bool physOnly) const
890 { 917 {
891 G4bool isFirst = true; 918 G4bool isFirst = true;
892 const G4ProcessVector* pv = fProcessManager- 919 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeGPIL);
893 G4int thisIdx(-1); 920 G4int thisIdx(-1);
894 for ( auto i = 0; i < (G4int)pv->size(); ++i << 921 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
895 { <<
896 if ( (*pv)(i) == this ) { thisIdx = i; bre <<
897 } <<
898 if ( thisIdx < 0 ) return false; // -- to ig 922 if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
899 for ( std::size_t i=0; i<(fSharedData->fBias << 923 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
900 { << 924 {
901 if ( (fSharedData->fBiasingProcessInterfac << 925 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
902 { << 926 {
903 G4int thatIdx(-1); << 927 G4int thatIdx(-1);
904 for (auto j = 0; j < (G4int)pv->size(); << 928 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
905 { << 929 if ( thatIdx >= 0 ) // -- to ignore pure along processes
906 if ( (*pv)(j) == (fSharedData->fBiasin << 930 {
907 { << 931 if ( thisIdx > thatIdx )
908 thatIdx = j; break; << 932 {
909 } << 933 isFirst = false;
910 } << 934 break;
911 if ( thatIdx >= 0 ) // -- to ignore pure << 935 }
912 { << 936 }
913 if ( thisIdx > thatIdx ) << 937 }
914 { <<
915 isFirst = false; <<
916 break; <<
917 } <<
918 } <<
919 } 938 }
920 } <<
921 return isFirst; 939 return isFirst;
922 } 940 }
923 941
924 G4bool G4BiasingProcessInterface:: << 942
925 IsLastPostStepGPILInterface(G4bool physOnly) c << 943 G4bool G4BiasingProcessInterface::IsLastPostStepGPILInterface(G4bool physOnly) const
926 { 944 {
927 G4bool isLast = true; 945 G4bool isLast = true;
928 const G4ProcessVector* pv = fProcessManager- 946 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeGPIL);
929 G4int thisIdx(-1); 947 G4int thisIdx(-1);
930 for (auto i = 0; i < (G4int)pv->size(); ++i << 948 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
931 { <<
932 if ( (*pv)(i) == this ) { thisIdx = i; bre <<
933 } <<
934 if ( thisIdx < 0 ) return false; // -- to ig 949 if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
935 for (std::size_t i=0; i<(fSharedData->fBiasi << 950 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
936 { << 951 {
937 if ( (fSharedData->fBiasingProcessInterfac << 952 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
938 { << 953 {
939 G4int thatIdx(-1); << 954 G4int thatIdx(-1);
940 for (auto j = 0; j < (G4int)pv->size(); << 955 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
941 { << 956 if ( thatIdx >= 0 ) // -- to ignore pure along processes
942 if ( (*pv)(j) == (fSharedData->fBiasin << 957 {
943 { << 958 if ( thisIdx < thatIdx )
944 thatIdx = j; break; << 959 {
945 } << 960 isLast = false;
946 } << 961 break;
947 if ( thatIdx >= 0 ) // -- to ignore pure << 962 }
948 { << 963 }
949 if ( thisIdx < thatIdx ) << 964 }
950 { <<
951 isLast = false; <<
952 break; <<
953 } <<
954 } <<
955 } 965 }
956 } <<
957 return isLast; 966 return isLast;
958 } 967 }
959 968
960 G4bool G4BiasingProcessInterface:: << 969
961 IsFirstPostStepDoItInterface(G4bool physOnly) << 970 G4bool G4BiasingProcessInterface::IsFirstPostStepDoItInterface(G4bool physOnly) const
962 { 971 {
963 G4bool isFirst = true; 972 G4bool isFirst = true;
964 const G4ProcessVector* pv = fProcessManager- 973 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeDoIt);
965 G4int thisIdx(-1); 974 G4int thisIdx(-1);
966 for (auto i = 0; i < (G4int)pv->size(); ++i << 975 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
967 { <<
968 if ( (*pv)(i) == this ) { thisIdx = i; bre <<
969 } <<
970 if ( thisIdx < 0 ) return false; // -- to ig 976 if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
971 for (std::size_t i=0; i<(fSharedData->fBiasi << 977 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
972 { << 978 {
973 if ( (fSharedData->fBiasingProcessInterfac << 979 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
974 { << 980 {
975 G4int thatIdx(-1); << 981 G4int thatIdx(-1);
976 for (auto j = 0; j < (G4int)pv->size(); << 982 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
977 { << 983 if ( thatIdx >= 0 ) // -- to ignore pure along processes
978 if ( (*pv)(j) == (fSharedData->fBiasin << 984 {
979 { << 985 if ( thisIdx > thatIdx )
980 thatIdx = j; break; << 986 {
981 } << 987 isFirst = false;
982 } << 988 break;
983 if ( thatIdx >= 0 ) // -- to ignore pure << 989 }
984 { << 990 }
985 if ( thisIdx > thatIdx ) << 991 }
986 { <<
987 isFirst = false; <<
988 break; <<
989 } <<
990 } <<
991 } 992 }
992 } <<
993 return isFirst; 993 return isFirst;
994 } 994 }
995 995
996 G4bool G4BiasingProcessInterface:: << 996
997 IsLastPostStepDoItInterface(G4bool physOnly) c << 997 G4bool G4BiasingProcessInterface::IsLastPostStepDoItInterface(G4bool physOnly) const
998 { 998 {
999 G4bool isLast = true; 999 G4bool isLast = true;
1000 const G4ProcessVector* pv = fProcessManager 1000 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeDoIt);
1001 G4int thisIdx(-1); 1001 G4int thisIdx(-1);
1002 for (auto i = 0; i < (G4int)pv->size(); ++i << 1002 for (G4int i = 0; i < pv->size(); i++ ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
1003 { <<
1004 if ( (*pv)(i) == this ) { thisIdx = i; br <<
1005 } <<
1006 if ( thisIdx < 0 ) return false; // -- to i 1003 if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
1007 for (std::size_t i=0; i<(fSharedData->fBias << 1004 for ( size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); i++ )
1008 { << 1005 {
1009 if ( (fSharedData->fBiasingProcessInterfa << 1006 if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
1010 { << 1007 {
1011 G4int thatIdx(-1); << 1008 G4int thatIdx(-1);
1012 for (auto j = 0; j < (G4int)pv->size(); << 1009 for (G4int j = 0; j < pv->size(); j++ ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
1013 { << 1010 if ( thatIdx >= 0 ) // -- to ignore pure along processes
1014 if ( (*pv)(j) == (fSharedData->fBiasi << 1011 {
1015 { << 1012 if ( thisIdx < thatIdx )
1016 thatIdx = j; break; << 1013 {
1017 } << 1014 isLast = false;
1018 } << 1015 break;
1019 if ( thatIdx >= 0 ) // -- to ignore pur << 1016 }
1020 { << 1017 }
1021 if ( thisIdx < thatIdx ) << 1018 }
1022 { <<
1023 isLast = false; <<
1024 break; <<
1025 } <<
1026 } <<
1027 } 1019 }
1028 } <<
1029 return isLast; 1020 return isLast;
1030 } 1021 }
1031 1022
>> 1023
1032 void G4BiasingProcessInterface::SetUpFirstLas 1024 void G4BiasingProcessInterface::SetUpFirstLastFlags()
1033 { 1025 {
1034 for (G4int iPhys = 0; iPhys < 2; ++iPhys) << 1026 for ( G4int iPhys = 0; iPhys < 2; iPhys++ )
1035 { << 1027 {
1036 G4bool physOnly = ( iPhys == 1 ); << 1028 G4bool physOnly = ( iPhys == 1 );
1037 fFirstLastFlags[IdxFirstLast( 1, 1, iPhys << 1029 fFirstLastFlags[IdxFirstLast( 1, 1, iPhys)] = IsFirstPostStepGPILInterface(physOnly);
1038 fFirstLastFlags[IdxFirstLast( 0, 1, iPhys << 1030 fFirstLastFlags[IdxFirstLast( 0, 1, iPhys)] = IsLastPostStepGPILInterface(physOnly);
1039 fFirstLastFlags[IdxFirstLast( 1, 0, iPhys << 1031 fFirstLastFlags[IdxFirstLast( 1, 0, iPhys)] = IsFirstPostStepDoItInterface(physOnly);
1040 fFirstLastFlags[IdxFirstLast( 0, 0, iPhys << 1032 fFirstLastFlags[IdxFirstLast( 0, 0, iPhys)] = IsLastPostStepDoItInterface(physOnly);
1041 } << 1033 }
1042 1034
1043 // -- for itself, for optimization: 1035 // -- for itself, for optimization:
1044 fIamFirstGPIL = GetIsFirstPostStepGPILInte 1036 fIamFirstGPIL = GetIsFirstPostStepGPILInterface( false );
1045 } 1037 }
1046 1038
>> 1039
1047 void G4BiasingProcessInterface::ResetForUnbia 1040 void G4BiasingProcessInterface::ResetForUnbiasedTracking()
1048 { 1041 {
1049 fOccurenceBiasingOperation = nullptr; << 1042 fOccurenceBiasingOperation = 0;
1050 fFinalStateBiasingOperation = nullptr; << 1043 fFinalStateBiasingOperation = 0;
1051 fNonPhysicsBiasingOperation = nullptr; << 1044 fNonPhysicsBiasingOperation = 0;
1052 fBiasingInteractionLaw = nullptr; << 1045 fBiasingInteractionLaw = 0;
1053 } 1046 }
1054 1047
1055 void G4BiasingProcessInterface:: << 1048
1056 InvokeWrappedProcessPostStepGPIL( const G4Tra << 1049 void G4BiasingProcessInterface::InvokeWrappedProcessPostStepGPIL( const G4Track& track,
1057 G4double pr << 1050 G4double previousStepSize,
1058 G4ForceCond << 1051 G4ForceCondition* condition )
1059 { 1052 {
1060 G4double usedPreviousStepSize = previousSte 1053 G4double usedPreviousStepSize = previousStepSize;
1061 // -- if the physics process has been under << 1054 // -- 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 1055 // -- 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 1056 // -- step. The pity is that PostStepGPIL and interaction length (cross-section)
1064 // -- calculations are done both in the Pos 1057 // -- calculations are done both in the PostStepGPIL of the process, while here we
1065 // -- are just interested in the calculatio 1058 // -- are just interested in the calculation of the cross-section. This is a pity
1066 // -- as this forces to re-generated a rand 1059 // -- as this forces to re-generated a random number for nothing.
1067 if ( fResetWrappedProcessInteractionLength 1060 if ( fResetWrappedProcessInteractionLength )
1068 { << 1061 {
1069 fResetWrappedProcessInteractionLength = f << 1062 fResetWrappedProcessInteractionLength = false;
1070 fWrappedProcess->ResetNumberOfInteraction << 1063 fWrappedProcess->ResetNumberOfInteractionLengthLeft();
1071 // -- We set "previous step size" as 0.0, << 1064 // -- We set "previous step size" as 0.0, to let the process believe this is first step:
1072 usedPreviousStepSize = 0.0; << 1065 usedPreviousStepSize = 0.0;
1073 } << 1066 }
1074 // -- GPIL response: 1067 // -- GPIL response:
1075 fWrappedProcessPostStepGPIL = fWrappedProce << 1068 fWrappedProcessPostStepGPIL = fWrappedProcess->PostStepGetPhysicalInteractionLength(track, usedPreviousStepSize, condition);
1076 fWrappedProcessForceCondition = *condition; << 1069 fWrappedProcessForceCondition = *condition;
1077 // -- and (inverse) cross-section: 1070 // -- and (inverse) cross-section:
1078 fWrappedProcessInteractionLength = fWrapped 1071 fWrappedProcessInteractionLength = fWrappedProcess->GetCurrentInteractionLength();
1079 } 1072 }
1080 1073
>> 1074
1081 void G4BiasingProcessInterface::ReorderBiasin 1075 void G4BiasingProcessInterface::ReorderBiasingVectorAsGPIL()
1082 { 1076 {
1083 // -- re-order vector of processes to match 1077 // -- re-order vector of processes to match that of the GPIL:
1084 std::vector < G4BiasingProcessInterface* > 1078 std::vector < G4BiasingProcessInterface* > tmpProcess ( fSharedData->fBiasingProcessInterfaces );
1085 ( fSharedData -> fBiasingProcessInterfaces 1079 ( fSharedData -> fBiasingProcessInterfaces ) . clear();
1086 ( fSharedData -> fPhysicsBiasingProcessInte 1080 ( fSharedData -> fPhysicsBiasingProcessInterfaces ) . clear();
1087 ( fSharedData -> fNonPhysicsBiasingProcessI 1081 ( fSharedData -> fNonPhysicsBiasingProcessInterfaces ) . clear();
1088 ( fSharedData -> fPublicBiasingProcessInter 1082 ( fSharedData -> fPublicBiasingProcessInterfaces ) . clear();
1089 ( fSharedData -> fPublicPhysicsBiasingProce 1083 ( fSharedData -> fPublicPhysicsBiasingProcessInterfaces ) . clear();
1090 ( fSharedData -> fPublicNonPhysicsBiasingPr 1084 ( fSharedData -> fPublicNonPhysicsBiasingProcessInterfaces ) . clear();
1091 1085
1092 const G4ProcessVector* pv = fProcessManager 1086 const G4ProcessVector* pv = fProcessManager->GetPostStepProcessVector(typeGPIL);
1093 for (auto i = 0; i < (G4int)pv->size(); ++i << 1087 for (G4int i = 0; i < pv->size(); i++ )
1094 { << 1088 {
1095 for (std::size_t j = 0; j < tmpProcess.si << 1089 for ( size_t j = 0; j < tmpProcess.size(); j++ )
1096 { << 1090 {
1097 if ( (*pv)(i) == tmpProcess[j] ) << 1091 if ( (*pv)(i) == tmpProcess[j] )
1098 { << 1092 {
1099 ( fSharedData->fBiasingProcessInterfa << 1093 ( fSharedData -> fBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1100 ( fSharedData->fPublicBiasingProcessI << 1094 ( fSharedData -> fPublicBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1101 if ( tmpProcess[j] -> fIsPhysicsBased << 1095 if ( tmpProcess[j] -> fIsPhysicsBasedBiasing )
1102 { << 1096 {
1103 ( fSharedData->fPhysicsBiasingProce << 1097 ( fSharedData -> fPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1104 ( fSharedData->fPublicPhysicsBiasin << 1098 ( fSharedData -> fPublicPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1105 } << 1099 }
1106 else << 1100 else
1107 { << 1101 {
1108 ( fSharedData -> fNonPhysicsBiasing << 1102 ( fSharedData -> fNonPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1109 ( fSharedData -> fPublicNonPhysicsB << 1103 ( fSharedData -> fPublicNonPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1110 } << 1104 }
1111 break; << 1105 break;
1112 } << 1106 }
>> 1107 }
1113 } 1108 }
1114 } <<
1115 } 1109 }
1116 1110