ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BOptrForceCollision.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BOptrForceCollision.cc
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 {
55  fCloningOperation = new G4BOptnCloning("Cloning");
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 {
79  fCloningOperation = new G4BOptnCloning("Cloning");
81 }
82 
83 
85 {
86  for ( std::map< const G4BiasingProcessInterface*, G4BOptnForceFreeFlight* >::iterator it = fFreeFlightOperations.begin() ;
87  it != fFreeFlightOperations.end() ;
88  it++ ) delete (*it).second;
90  delete fCloningOperation;
91 }
92 
93 
95 {
96  // -- Create ID for force collision:
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 
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):
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:
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:
191  else
192  {
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. ]
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. ]
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.
215  {
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):
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 
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:
272  if ( fCurrentTrackData != nullptr )
273  {
275  {
276  // -- takes "ownership" (some track data created before, left free, reuse it):
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  {
288  }
289  fCurrentTrackData->fForceCollisionState = ForceCollisionState::toBeCloned;
290  fInitialTrackWeight = track->GetWeight();
292  return fCloningOperation;
293  }
294 
295  // --
296  return nullptr;
297 }
298 
299 
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 {
312  fCurrentTrackData = nullptr;
313 }
314 
315 
317 {
318  // -- check for consistency, operator should have cleaned the track:
319  if ( fCurrentTrackData != nullptr )
320  {
322  {
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;
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  }
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