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