Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BOptrForceCollision.cc
Go to the documentation of this file.
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 #include "G4BOptrForceCollision.hh"
29 #include "G4PhysicsModelCatalog.hh"
30 
34 #include "G4BOptnCloning.hh"
35 
36 #include "G4Step.hh"
37 #include "G4StepPoint.hh"
38 #include "G4VProcess.hh"
39 
40 #include "G4ParticleDefinition.hh"
41 #include "G4ParticleTable.hh"
42 
43 #include "G4SystemOfUnits.hh"
44 
45 // -- §§ consider calling other constructor, thanks to C++11
47  : G4VBiasingOperator(name),
48  fForceCollisionModelID(-1),
49  fCurrentTrack(nullptr),
50  fCurrentTrackData(nullptr),
51  fInitialTrackWeight(-1.0),
52  fSetup(true)
53 {
54  fSharedForceInteractionOperation = new G4BOptnForceCommonTruncatedExp("SharedForceInteraction");
55  fCloningOperation = new G4BOptnCloning("Cloning");
56  fParticleToBias = G4ParticleTable::GetParticleTable()->FindParticle(particleName);
57 
58  if ( fParticleToBias == 0 )
59  {
61  ed << " Particle `" << particleName << "' not found !" << G4endl;
62  G4Exception(" G4BOptrForceCollision::G4BOptrForceCollision(...)",
63  "BIAS.GEN.07",
65  ed);
66  }
67 }
68 
69 
71  : G4VBiasingOperator(name),
72  fForceCollisionModelID(-1),
73  fCurrentTrack(nullptr),
74  fCurrentTrackData(nullptr),
75  fInitialTrackWeight(-1.0),
76  fSetup(true)
77 {
78  fSharedForceInteractionOperation = new G4BOptnForceCommonTruncatedExp("SharedForceInteraction");
79  fCloningOperation = new G4BOptnCloning("Cloning");
80  fParticleToBias = particle;
81 }
82 
83 
85 {
86  for ( std::map< const G4BiasingProcessInterface*, G4BOptnForceFreeFlight* >::iterator it = fFreeFlightOperations.begin() ;
87  it != fFreeFlightOperations.end() ;
88  it++ ) delete (*it).second;
89  delete fSharedForceInteractionOperation;
90  delete fCloningOperation;
91 }
92 
93 
95 {
96  // -- Create ID for force collision:
97  fForceCollisionModelID = G4PhysicsModelCatalog::Register("GenBiasForceCollision");
98 
99  // -- build free flight operations:
101 }
102 
103 
105 {
106  // -- start by remembering processes under biasing, and create needed biasing operations:
107  if ( fSetup )
108  {
109  // -- get back ID created in master thread in case of MT (or reget just created ID in sequential mode):
110  fForceCollisionModelID = G4PhysicsModelCatalog::Register("GenBiasForceCollision");
111 
112  const G4ProcessManager* processManager = fParticleToBias->GetProcessManager();
113  const G4BiasingProcessSharedData* interfaceProcessSharedData = G4BiasingProcessInterface::GetSharedData( processManager );
114  if ( interfaceProcessSharedData ) // -- sharedData tested, as is can happen a user attaches an operator
115  // -- to a volume but without defining BiasingProcessInterface processes.
116  {
117  for ( size_t i = 0 ; i < (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces()).size(); i++ )
118  {
119  const G4BiasingProcessInterface* wrapperProcess =
120  (interfaceProcessSharedData->GetPhysicsBiasingProcessInterfaces())[i];
121  G4String operationName = "FreeFlight-"+wrapperProcess->GetWrappedProcess()->GetProcessName();
122  fFreeFlightOperations[wrapperProcess] = new G4BOptnForceFreeFlight(operationName);
123  }
124  }
125  fSetup = false;
126  }
127 }
128 
129 
131 {
132 }
133 
134 
135 G4VBiasingOperation* G4BOptrForceCollision::ProposeOccurenceBiasingOperation(const G4Track* track, const G4BiasingProcessInterface* callingProcess)
136 {
137  // -- does nothing if particle is not of requested type:
138  if ( track->GetDefinition() != fParticleToBias ) return 0;
139 
140  // -- trying to get auxiliary track data...
141  if ( fCurrentTrackData == nullptr )
142  {
143  // ... and if the track has no aux. track data, it means its biasing is not started yet (note that cloning is to happen first):
144  fCurrentTrackData = (G4BOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID));
145  if ( fCurrentTrackData == nullptr ) return nullptr;
146  }
147 
148 
149  // -- Send force free flight to the callingProcess:
150  // ------------------------------------------------
151  // -- The track has been cloned in the previous step, it has now to be
152  // -- forced for a free flight.
153  // -- This track will fly with 0.0 weight during its forced flight:
154  // -- it is to forbid the double counting with the force interaction track.
155  // -- Its weight is restored at the end of its free flight, this weight
156  // -- being its initial weight * the weight for the free flight travel,
157  // -- this last one being per process. The initial weight is common, and is
158  // -- arbitrary asked to the first operation to take care of it.
159  if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight )
160  {
161  G4BOptnForceFreeFlight* operation = fFreeFlightOperations[callingProcess];
162  if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. )
163  {
164  // -- the initial track weight will be restored only by the first DoIt free flight:
165  operation->ResetInitialTrackWeight(fInitialTrackWeight);
166  return operation;
167  }
168  else
169  {
170  return nullptr;
171  }
172  }
173 
174 
175  // -- Send force interaction operation to the callingProcess:
176  // ----------------------------------------------------------
177  // -- at this level, a copy of the track entering the volume was
178  // -- generated (borned) earlier. This copy will make the forced
179  // -- interaction in the volume.
180  if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )
181  {
182  // -- Remember if this calling process is the first of the physics wrapper in
183  // -- the PostStepGPIL loop (using default argument of method below):
184  G4bool isFirstPhysGPIL = callingProcess-> GetIsFirstPostStepGPILInterface();
185 
186  // -- [*first process*] Initialize or update force interaction operation:
187  if ( isFirstPhysGPIL )
188  {
189  // -- first step of cloned track, initialize the forced interaction operation:
190  if ( track->GetCurrentStepNumber() == 1 ) fSharedForceInteractionOperation->Initialize( track );
191  else
192  {
193  if ( fSharedForceInteractionOperation->GetInitialMomentum() != track->GetMomentum() )
194  {
195  // -- means that some other physics process, not under control of the forced interaction operation,
196  // -- has occured, need to re-initialize the operation as distance to boundary has changed.
197  // -- [ Note the re-initialization is only possible for a Markovian law. ]
198  fSharedForceInteractionOperation->Initialize( track );
199  }
200  else
201  {
202  // -- means that some other non-physics process (biasing or not, like step limit), has occured,
203  // -- but track conserves its momentum direction, only need to reduced the maximum distance for
204  // -- forced interaction.
205  // -- [ Note the update is only possible for a Markovian law. ]
206  fSharedForceInteractionOperation->UpdateForStep( track->GetStep() );
207  }
208  }
209  }
210 
211  // -- [*all processes*] Sanity check : it may happen in limit cases that distance to
212  // -- out is zero, weight would be infinite in this case: abort forced interaction
213  // -- and abandon biasing.
214  if ( fSharedForceInteractionOperation->GetMaximumDistance() < DBL_MIN )
215  {
216  fCurrentTrackData->Reset();
217  return 0;
218  }
219 
220  // -- [* first process*] collect cross-sections and sample force law to determine interaction length
221  // -- and winning process:
222  if ( isFirstPhysGPIL )
223  {
224  // -- collect cross-sections:
225  // -- ( Remember that the first of the G4BiasingProcessInterface triggers the update
226  // -- of these cross-sections )
227  const G4BiasingProcessSharedData* sharedData = callingProcess->GetSharedData();
228  for ( size_t i = 0 ; i < (sharedData->GetPhysicsBiasingProcessInterfaces()).size() ; i++ )
229  {
230  const G4BiasingProcessInterface* wrapperProcess = ( sharedData->GetPhysicsBiasingProcessInterfaces() )[i];
231  G4double interactionLength = wrapperProcess->GetWrappedProcess()->GetCurrentInteractionLength();
232  // -- keep only well defined cross-sections, other processes are ignored. These are not pathological cases
233  // -- but cases where a threhold effect par example (eg pair creation) may be at play. (**!**)
234  if ( interactionLength < DBL_MAX/10. )
235  fSharedForceInteractionOperation->AddCrossSection( wrapperProcess->GetWrappedProcess(), 1.0/interactionLength );
236  }
237  // -- sample the shared law (interaction length, and winning process):
238  if ( fSharedForceInteractionOperation->GetNumberOfSharing() > 0 ) fSharedForceInteractionOperation->Sample();
239  }
240 
241  // -- [*all processes*] Send operation for processes with well defined XS (see "**!**" above):
242  G4VBiasingOperation* operationToReturn = nullptr;
243  if ( callingProcess->GetWrappedProcess()->GetCurrentInteractionLength() < DBL_MAX/10. ) operationToReturn = fSharedForceInteractionOperation;
244  return operationToReturn;
245 
246 
247  } // -- end of "if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )"
248 
249 
250  // -- other cases here: particle appearing in the volume by some
251  // -- previous interaction : we decide to not bias these.
252  return 0;
253 
254 }
255 
256 
257 G4VBiasingOperation* G4BOptrForceCollision::ProposeNonPhysicsBiasingOperation(const G4Track* track,
258  const G4BiasingProcessInterface* /* callingProcess */)
259 {
260  if ( track->GetDefinition() != fParticleToBias ) return nullptr;
261 
262  // -- Apply biasing scheme only to tracks entering the volume.
263  // -- A "cloning" is done:
264  // -- - the primary will be forced flight under a zero weight up to volume exit,
265  // -- where the weight will be restored with proper weight for free flight
266  // -- - the clone will be forced to interact in the volume.
267  if ( track->GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary ) // -- §§§ extent to case of a track shoot on the boundary ?
268  {
269  // -- check that track is free of undergoing biasing scheme ( no biasing data, or no active active )
270  // -- Get possible track data:
271  fCurrentTrackData = (G4BOptrForceCollisionTrackData*)(track->GetAuxiliaryTrackInformation(fForceCollisionModelID));
272  if ( fCurrentTrackData != nullptr )
273  {
274  if ( fCurrentTrackData->IsFreeFromBiasing() )
275  {
276  // -- takes "ownership" (some track data created before, left free, reuse it):
277  fCurrentTrackData->fForceCollisionOperator = this ;
278  }
279  else
280  {
281  // §§§ Would something be really wrong in this case ? Could this be that a process made a zero step ?
282  }
283  }
284  else
285  {
286  fCurrentTrackData = new G4BOptrForceCollisionTrackData( this );
287  track->SetAuxiliaryTrackInformation(fForceCollisionModelID, fCurrentTrackData);
288  }
289  fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeCloned;
290  fInitialTrackWeight = track->GetWeight();
291  fCloningOperation->SetCloneWeights(0.0, fInitialTrackWeight);
292  return fCloningOperation;
293  }
294 
295  // --
296  return nullptr;
297 }
298 
299 
300 G4VBiasingOperation* G4BOptrForceCollision::ProposeFinalStateBiasingOperation(const G4Track*, const G4BiasingProcessInterface* callingProcess)
301 {
302  // -- Propose at final state generation the same operation which was proposed at GPIL level,
303  // -- (which is either the force free flight or the force interaction operation).
304  // -- count on the process interface to collect this operation:
305  return callingProcess->GetCurrentOccurenceBiasingOperation();
306 }
307 
308 
310 {
311  fCurrentTrack = track;
312  fCurrentTrackData = nullptr;
313 }
314 
315 
317 {
318  // -- check for consistency, operator should have cleaned the track:
319  if ( fCurrentTrackData != nullptr )
320  {
321  if ( !fCurrentTrackData->IsFreeFromBiasing() )
322  {
323  if ( (fCurrentTrack->GetTrackStatus() == fStopAndKill) || (fCurrentTrack->GetTrackStatus() == fKillTrackAndSecondaries) )
324  {
326  ed << "Current track deleted while under biasing by " << GetName() << ". Will result in inconsistencies.";
327  G4Exception(" G4BOptrForceCollision::EndTracking()",
328  "BIAS.GEN.18",
329  JustWarning,
330  ed);
331  }
332  }
333  }
334 }
335 
336 
339  G4VBiasingOperation* operationApplied,
340  const G4VParticleChange* )
341 {
342 
343  if ( fCurrentTrackData == nullptr )
344  {
345  if ( BAC != BAC_None )
346  {
348  ed << " Internal inconsistency : please submit bug report. " << G4endl;
349  G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
350  "BIAS.GEN.20.1",
351  JustWarning,
352  ed);
353  }
354  return;
355  }
356 
357  if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeCloned )
358  {
359  fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeFreeFlight;
360  auto cloneData = new G4BOptrForceCollisionTrackData( this );
361  cloneData->fForceCollisionState = ForceCollisionState::toBeForced;
362  fCloningOperation->GetCloneTrack()->SetAuxiliaryTrackInformation(fForceCollisionModelID, cloneData);
363  }
364  else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeFreeFlight )
365  {
366  if ( fFreeFlightOperations[callingProcess]->OperationComplete() ) fCurrentTrackData->Reset(); // -- off biasing for this track
367  }
368  else if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )
369  {
370  if ( operationApplied != fSharedForceInteractionOperation )
371  {
373  ed << " Internal inconsistency : please submit bug report. " << G4endl;
374  G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
375  "BIAS.GEN.20.2",
376  JustWarning,
377  ed);
378  }
379  if ( fSharedForceInteractionOperation->GetInteractionOccured() )
380  {
381  if ( operationApplied != fSharedForceInteractionOperation )
382  {
384  ed << " Internal inconsistency : please submit bug report. " << G4endl;
385  G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
386  "BIAS.GEN.20.3",
387  JustWarning,
388  ed);
389  }
390  }
391  }
392  else
393  {
394  if ( fCurrentTrackData->fForceCollisionState != ForceCollisionState::free )
395  {
397  ed << " Internal inconsistency : please submit bug report. " << G4endl;
398  G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
399  "BIAS.GEN.20.4",
400  JustWarning,
401  ed);
402  }
403  }
404 }
405 
406 
408  G4VBiasingOperation* /*occurenceOperationApplied*/, G4double /*weightForOccurenceInteraction*/,
409  G4VBiasingOperation* finalStateOperationApplied, const G4VParticleChange* /*particleChangeProduced*/ )
410 {
411 
412  if ( fCurrentTrackData->fForceCollisionState == ForceCollisionState::toBeForced )
413  {
414  if ( finalStateOperationApplied != fSharedForceInteractionOperation )
415  {
417  ed << " Internal inconsistency : please submit bug report. " << G4endl;
418  G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
419  "BIAS.GEN.20.5",
420  JustWarning,
421  ed);
422  }
423  if ( fSharedForceInteractionOperation->GetInteractionOccured() ) fCurrentTrackData->Reset(); // -- off biasing for this track
424  }
425  else
426  {
428  ed << " Internal inconsistency : please submit bug report. " << G4endl;
429  G4Exception(" G4BOptrForceCollision::OperationApplied(...)",
430  "BIAS.GEN.20.6",
431  JustWarning,
432  ed);
433  }
434 
435 }
436 
G4ParticleDefinition * GetDefinition() const
const XML_Char * name
Definition: expat.h:151
G4VBiasingOperation * GetCurrentOccurenceBiasingOperation() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4Track * GetCloneTrack() const
void SetAuxiliaryTrackInformation(G4int idx, G4VAuxiliaryTrackInformation *info) const
Definition: G4Track.cc:324
G4StepStatus GetStepStatus() const
void AddCrossSection(const G4VProcess *, G4double)
G4TrackStatus GetTrackStatus() const
const std::vector< const G4BiasingProcessInterface * > & GetPhysicsBiasingProcessInterfaces() const
virtual void EndTracking() final
const G4ThreeVector & GetInitialMomentum() const
const G4Step * GetStep() const
G4BOptrForceCollision(G4String particleToForce, G4String name="ForceCollision")
G4BiasingAppliedCase
G4VProcess * GetWrappedProcess() const
void OperationApplied(const G4BiasingProcessInterface *callingProcess, G4BiasingAppliedCase biasingCase, G4VBiasingOperation *operationApplied, const G4VParticleChange *particleChangeProduced) final
G4StepPoint * GetPreStepPoint() const
void ResetInitialTrackWeight(G4double w)
G4int GetCurrentStepNumber() const
virtual void StartRun() final
bool G4bool
Definition: G4Types.hh:79
G4double GetCurrentInteractionLength() const
Definition: G4VProcess.hh:462
virtual void StartTracking(const G4Track *track) final
virtual void Configure() final
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4VAuxiliaryTrackInformation * GetAuxiliaryTrackInformation(G4int idx) const
Definition: G4Track.cc:340
void SetCloneWeights(G4double clone1Weight, G4double clone2Weight)
G4ThreeVector GetMomentum() const
static G4ParticleTable * GetParticleTable()
G4ProcessManager * GetProcessManager() const
static G4int Register(const G4String &)
#define DBL_MIN
Definition: templates.hh:75
G4double GetWeight() const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83
virtual void ConfigureForWorker() final
const G4String GetName() const
const G4BiasingProcessSharedData * GetSharedData() const