| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4BOptrForceCollision
27 // --------------------------------------------------------------------
28
29 #include "G4BOptrForceCollision.hh"
30 #include "G4BOptrForceCollisionTrackData.hh"
31 #include "G4BiasingProcessInterface.hh"
32 #include "G4PhysicsModelCatalog.hh"
33
34 #include "G4BOptnForceCommonTruncatedExp.hh"
35 #include "G4ILawCommonTruncatedExp.hh"
36 #include "G4BOptnForceFreeFlight.hh"
37 #include "G4BOptnCloning.hh"
38
39 #include "G4Step.hh"
40 #include "G4StepPoint.hh"
41 #include "G4VProcess.hh"
42
43 #include "G4ParticleDefinition.hh"
44 #include "G4ParticleTable.hh"
45
46 #include "G4SystemOfUnits.hh"
47
48 // -- Consider calling other constructor, thanks to C++11
49 G4BOptrForceCollision::
50 G4BOptrForceCollision(const G4String& particleName,
51 const G4String& name)
52 : G4VBiasingOperator(name),
53 fForceCollisionModelID(G4PhysicsModelCatalog::GetModelID("model_GenBiasForceCollision"))
54 {
55 fSharedForceInteractionOperation = new G4BOptnForceCommonTruncatedExp("SharedForceInteraction");
56 fCloningOperation = new G4BOptnCloning("Cloning");
57 fParticleToBias = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
58
59 if ( fParticleToBias == nullptr )
60 {
61 G4ExceptionDescription ed;
62 ed << " Particle `" << particleName << "' not found !" << G4endl;
63 G4Exception(" G4BOptrForceCollision::G4BOptrForceCollision(...)",
64 "BIAS.GEN.07", JustWarning, ed);
65 }
66 }
67
68 G4BOptrForceCollision::
69 G4BOptrForceCollision(const G4ParticleDefinition* particle,
70 const G4String& name)
71 : G4VBiasingOperator(name),
72 fForceCollisionModelID(G4PhysicsModelCatalog::GetModelID("model_GenBiasForceCollision"))
73 {
74 fSharedForceInteractionOperation = new G4BOptnForceCommonTruncatedExp("SharedForceInteraction");
75 fCloningOperation = new G4BOptnCloning("Cloning");
76 fParticleToBias = particle;
77 }
78
79 G4BOptrForceCollision::~G4BOptrForceCollision()
80 {
81 for ( auto it = fFreeFlightOperations.cbegin();
82 it != fFreeFlightOperations.cend(); ++it )
83 {
84 delete (*it).second;
85 }
86 delete fSharedForceInteractionOperation;
87 delete fCloningOperation;
88 }
89
90 void G4BOptrForceCollision::Configure()
91 {
92 // -- build free flight operations:
93 ConfigureForWorker();
94 }
95
96 void G4BOptrForceCollision::ConfigureForWorker()
97 {
98 // -- start by remembering processes under biasing,
99 // and create needed biasing operations:
100 if ( fSetup )
101 {
102 const G4ProcessManager* processManager = fParticleToBias->GetProcessManager();
103 const G4BiasingProcessSharedData* interfaceProcessSharedData = G4BiasingProcessInterface::GetSharedData( processManager );
104 if ( interfaceProcessSharedData ) // -- sharedData tested, as is can happen a user attaches an operator
105 // -- to a volume but without defining BiasingProcessInterface processes.
106 {
107 for ( std::size_t i = 0 ; i < (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces()).size(); ++i )
108 {
109 const G4BiasingProcessInterface* wrapperProcess =
110 (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces())[i];
111 const G4String& operationName = "FreeFlight-"+wrapperProcess->GetWrappedProcess()->GetProcessName();
112 fFreeFlightOperations[wrapperProcess] = new G4BOptnForceFreeFlight(operationName);
113 }
114 }
115 fSetup = false;
116 }
117 }
118
119 void G4BOptrForceCollision::StartRun()
120 {
121 }
122
123 G4VBiasingOperation*
124 G4BOptrForceCollision::ProposeOccurenceBiasingOperation(const G4Track* track,
125 const G4BiasingProcessInterface* callingProcess)
126 {
127 // -- does nothing if particle is not of requested type:
128 if ( track->GetDefinition() != fParticleToBias ) return nullptr;
129
130 // -- trying to get auxiliary track data...
131 if ( fCurrentTrackData == nullptr )
132 {
133 // ... and if the track has no aux. track data, it means its biasing is not started yet (note that cloning is to happen first):
134 fCurrentTrackData = (G4BOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID));
135 if ( fCurrentTrackData == nullptr ) return nullptr;
136 }
137
138 // -- Send force free flight to the callingProcess:
139 // ------------------------------------------------
140 // -- The track has been cloned in the previous step, it has now to be
141 // -- forced for a free flight.
142 // -- This track will fly with 0.0 weight during its forced flight:
143 // -- it is to forbid the double counting with the force interaction track.
144 // -- Its weight is restored at the end of its free flight, this weight
145 // -- being its initial weight * the weight for the free flight travel,
146 // -- this last one being per process. The initial weight is common, and is
147 // -- arbitrary asked to the first operation to take care of it.
148 if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight )
149 {
150 G4BOptnForceFreeFlight* operation = fFreeFlightOperations[callingProcess];
151 if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. )
152 {
153 // -- the initial track weight will be restored only by the first DoIt free flight:
154 operation->ResetInitialTrackWeight(fInitialTrackWeight);
155 return operation;
156 }
157 else
158 {
159 return nullptr;
160 }
161 }
162
163 // -- Send force interaction operation to the callingProcess:
164 // ----------------------------------------------------------
165 // -- at this level, a copy of the track entering the volume was
166 // -- generated (borned) earlier. This copy will make the forced
167 // -- interaction in the volume.
168 if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )
169 {
170 // -- Remember if this calling process is the first of the physics wrapper in
171 // -- the PostStepGPIL loop (using default argument of method below):
172 G4bool isFirstPhysGPIL = callingProcess-> GetIsFirstPostStepGPILInterface();
173
174 // -- [*first process*] Initialize or update force interaction operation:
175 if ( isFirstPhysGPIL )
176 {
177 // -- first step of cloned track, initialize the forced interaction operation:
178 if ( track->GetCurrentStepNumber() == 1 )
179 {
180 fSharedForceInteractionOperation->Initialize( track );
181 }
182 else
183 {
184 if ( fSharedForceInteractionOperation->GetInitialMomentum() != track->GetMomentum() )
185 {
186 // -- means that some other physics process, not under control of the forced interaction operation,
187 // -- has occured, need to re-initialize the operation as distance to boundary has changed.
188 // -- [ Note the re-initialization is only possible for a Markovian law. ]
189 fSharedForceInteractionOperation->Initialize( track );
190 }
191 else
192 {
193 // -- means that some other non-physics process (biasing or not, like step limit), has occured,
194 // -- but track conserves its momentum direction, only need to reduced the maximum distance for
195 // -- forced interaction.
196 // -- [ Note the update is only possible for a Markovian law. ]
197 fSharedForceInteractionOperation->UpdateForStep( track->GetStep() );
198 }
199 }
200 }
201
202 // -- [*all processes*] Sanity check : it may happen in limit cases that distance to
203 // -- out is zero, weight would be infinite in this case: abort forced interaction
204 // -- and abandon biasing.
205 if ( fSharedForceInteractionOperation->GetMaximumDistance() < DBL_MIN )
206 {
207 fCurrentTrackData->Reset();
208 return nullptr;
209 }
210
211 // -- [* first process*] collect cross-sections and sample force law to determine interaction length
212 // -- and winning process:
213 if ( isFirstPhysGPIL )
214 {
215 // -- collect cross-sections:
216 // -- ( Remember that the first of the G4BiasingProcessInterface triggers the update
217 // -- of these cross-sections )
218 const G4BiasingProcessSharedData* sharedData = callingProcess->GetSharedData();
219 for ( std::size_t i = 0 ; i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size() ; ++i )
220 {
221 const G4BiasingProcessInterface* wrapperProcess = ( sharedData->GetPhysicsBiasingProcessInterfaces() )[i];
222 G4double interactionLength = wrapperProcess->GetWrappedProcess()->GetCurrentInteractionLength();
223 // -- keep only well defined cross-sections, other processes are ignored. These are not pathological cases
224 // -- but cases where a threhold effect par example (eg pair creation) may be at play. (**!**)
225 if ( interactionLength < DBL_MAX/10. )
226 fSharedForceInteractionOperation->AddCrossSection( wrapperProcess->GetWrappedProcess(), 1.0/interactionLength );
227 }
228 // -- sample the shared law (interaction length, and winning process):
229 if ( fSharedForceInteractionOperation->GetNumberOfSharing() > 0 )
230 fSharedForceInteractionOperation->Sample();
231 }
232
233 // -- [*all processes*] Send operation for processes with well defined XS (see "**!**" above):
234 G4VBiasingOperation* operationToReturn = nullptr;
235 if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. )
236 operationToReturn = fSharedForceInteractionOperation;
237 return operationToReturn;
238
239 } // -- end of "if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )"
240
241
242 // -- other cases here: particle appearing in the volume by some
243 // -- previous interaction : we decide to not bias these.
244 return nullptr;
245 }
246
247 G4VBiasingOperation* G4BOptrForceCollision::
248 ProposeNonPhysicsBiasingOperation(const G4Track* track,
249 const G4BiasingProcessInterface* /* callingProcess */)
250 {
251 if ( track->GetDefinition() != fParticleToBias ) return nullptr;
252
253 // -- Apply biasing scheme only to tracks entering the volume.
254 // -- A "cloning" is done:
255 // -- - the primary will be forced flight under a zero weight up to volume exit,
256 // -- where the weight will be restored with proper weight for free flight
257 // -- - the clone will be forced to interact in the volume.
258 if ( track->GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ) // -- ยงยงยง extent to case of a track shoot on the boundary ?
259 {
260 // -- check that track is free of undergoing biasing scheme ( no biasing data, or no active active )
261 // -- Get possible track data:
262 fCurrentTrackData = (G4BOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID));
263 if ( fCurrentTrackData != nullptr )
264 {
265 if ( fCurrentTrackData->IsFreeFromBiasing() )
266 {
267 // -- takes "ownership" (some track data created before, left free, reuse it):
268 fCurrentTrackData->fForceCollisionOperator = this ;
269 }
270 else
271 {
272 // Would something be really wrong in this case ?
273 // Could this be that a process made a zero step ?
274 }
275 }
276 else
277 {
278 fCurrentTrackData = new G4BOptrForceCollisionTrackData( this );
279 track->SetAuxiliaryTrackInformation(fForceCollisionModelID, fCurrentTrackData);
280 }
281 fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeCloned;
282 fInitialTrackWeight = track->GetWeight();
283 fCloningOperation->SetCloneWeights(0.0, fInitialTrackWeight);
284
285 return fCloningOperation;
286 }
287
288 return nullptr;
289 }
290
291 G4VBiasingOperation* G4BOptrForceCollision::
292 ProposeFinalStateBiasingOperation(const G4Track*,
293 const G4BiasingProcessInterface* callingProcess)
294 {
295 // -- Propose at final state generation the same operation which was proposed at GPIL level,
296 // -- (which is either the force free flight or the force interaction operation).
297 // -- count on the process interface to collect this operation:
298 return callingProcess->GetCurrentOccurenceBiasingOperation();
299 }
300
301 void G4BOptrForceCollision::StartTracking( const G4Track* track )
302 {
303 fCurrentTrack = track;
304 fCurrentTrackData = nullptr;
305 }
306
307 void G4BOptrForceCollision::EndTracking()
308 {
309 // -- check for consistency, operator should have cleaned the track:
310 if ( fCurrentTrackData != nullptr )
311 {
312 if ( !fCurrentTrackData->IsFreeFromBiasing() )
313 {
314 if ( (fCurrentTrack->GetTrackStatus() == fStopAndKill)
315 || (fCurrentTrack->GetTrackStatus() == fKillTrackAndSecondaries) )
316 {
317 G4ExceptionDescription ed;
318 ed << "Current track deleted while under biasing by "
319 << GetName() << ". Will result in inconsistencies.";
320 G4Exception(" G4BOptrForceCollision::EndTracking()",
321 "BIAS.GEN.18", JustWarning, ed);
322 }
323 }
324 }
325 }
326
327 void G4BOptrForceCollision::
328 OperationApplied( const G4BiasingProcessInterface* callingProcess,
329 G4BiasingAppliedCase BAC,
330 G4VBiasingOperation* operationApplied,
331 const G4VParticleChange* )
332 {
333 if ( fCurrentTrackData == nullptr )
334 {
335 if ( BAC != BAC_None )
336 {
337 G4ExceptionDescription ed;
338 ed << " Internal inconsistency : please submit bug report. " << G4endl;
339 G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
340 "BIAS.GEN.20.1", JustWarning, ed);
341 }
342 return;
343 }
344
345 if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeCloned )
346 {
347 fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeFreeFlight;
348 auto cloneData = new G4BOptrForceCollisionTrackData( this );
349 cloneData->fForceCollisionState = ForceCollisionState::toBeForced;
350 fCloningOperation->GetCloneTrack()->SetAuxiliaryTrackInformation(fForceCollisionModelID, cloneData);
351 }
352 else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight )
353 {
354 if ( fFreeFlightOperations[callingProcess]->OperationComplete() )
355 fCurrentTrackData->Reset(); // -- off biasing for this track
356 }
357 else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )
358 {
359 if ( operationApplied != fSharedForceInteractionOperation )
360 {
361 G4ExceptionDescription ed;
362 ed << " Internal inconsistency : please submit bug report. " << G4endl;
363 G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
364 "BIAS.GEN.20.2", JustWarning, ed);
365 }
366 if ( fSharedForceInteractionOperation->GetInteractionOccured() )
367 {
368 if ( operationApplied != fSharedForceInteractionOperation )
369 {
370 G4ExceptionDescription ed;
371 ed << " Internal inconsistency : please submit bug report. " << G4endl;
372 G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
373 "BIAS.GEN.20.3", JustWarning, ed);
374 }
375 }
376 }
377 else
378 {
379 if ( fCurrentTrackData->fForceCollisionState != ForceCollisionState::free )
380 {
381 G4ExceptionDescription ed;
382 ed << " Internal inconsistency : please submit bug report. " << G4endl;
383 G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
384 "BIAS.GEN.20.4", JustWarning, ed);
385 }
386 }
387 }
388
389 void G4BOptrForceCollision::
390 OperationApplied( const G4BiasingProcessInterface* /*callingProcess*/,
391 G4BiasingAppliedCase /*biasingCase*/,
392 G4VBiasingOperation* /*occurenceOperationApplied*/,
393 G4double /*weightForOccurenceInteraction*/,
394 G4VBiasingOperation* finalStateOperationApplied,
395 const G4VParticleChange* /*particleChangeProduced*/ )
396 {
397 if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )
398 {
399 if ( finalStateOperationApplied != fSharedForceInteractionOperation )
400 {
401 G4ExceptionDescription ed;
402 ed << " Internal inconsistency : please submit bug report. " << G4endl;
403 G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
404 "BIAS.GEN.20.5", JustWarning, ed);
405 }
406 if ( fSharedForceInteractionOperation->GetInteractionOccured() )
407 fCurrentTrackData->Reset(); // -- off biasing for this track
408 }
409 else
410 {
411 G4ExceptionDescription ed;
412 ed << " Internal inconsistency : please submit bug report. " << G4endl;
413 G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
414 "BIAS.GEN.20.6", JustWarning, ed);
415 }
416 }
417