ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BiasingProcessInterface.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BiasingProcessInterface.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 //
27 #include "G4VBiasingOperator.hh"
28 #include "G4VBiasingOperation.hh"
33 #include "G4ProcessManager.hh"
34 #include "G4BiasingAppliedCase.hh"
36 
41 
43  : G4VProcess ( name ),
44  fCurrentTrack ( nullptr ),
45  fPreviousStepSize (-1.0), fCurrentMinimumStep( -1.0 ), fProposedSafety ( -1.0),
46  fOccurenceBiasingOperation( nullptr ), fFinalStateBiasingOperation( nullptr ), fNonPhysicsBiasingOperation( nullptr ),
47  fPreviousOccurenceBiasingOperation( nullptr ), fPreviousFinalStateBiasingOperation( nullptr ), fPreviousNonPhysicsBiasingOperation( nullptr ),
48  fResetWrappedProcessInteractionLength( true ),
49  fWrappedProcess ( nullptr ),
50  fIsPhysicsBasedBiasing ( false ),
51  fWrappedProcessIsAtRest ( false ),
52  fWrappedProcessIsAlong ( false ),
53  fWrappedProcessIsPost ( false ),
54  fWrappedProcessPostStepGPIL ( -1.0 ),
55  fBiasingPostStepGPIL ( -1.0 ),
56  fWrappedProcessInteractionLength ( -1.0 ),
57  fWrappedProcessForceCondition ( NotForced ),
58  fBiasingForceCondition ( NotForced ),
59  fWrappedProcessAlongStepGPIL ( -1.0 ),
60  fBiasingAlongStepGPIL ( -1.0 ),
61  fWrappedProcessGPILSelection ( NotCandidateForSelection ),
62  fBiasingGPILSelection ( NotCandidateForSelection ),
63  fBiasingInteractionLaw ( nullptr ),
64  fPreviousBiasingInteractionLaw ( nullptr ),
65  fPhysicalInteractionLaw ( nullptr ),
66  fOccurenceBiasingParticleChange ( nullptr ),
67  fDummyParticleChange ( nullptr ),
68  fIamFirstGPIL ( false ),
69  fProcessManager ( nullptr ),
70  fSharedData ( nullptr )
71 {
72  for (G4int i = 0 ; i < 8 ; i++) fFirstLastFlags[i] = false;
73  fResetInteractionLaws.Put( true );
74  fCommonStart .Put( true );
75  fCommonEnd .Put( true );
76  fDoCommonConfigure .Put( true );
77 }
78 
79 
81  G4bool wrappedIsAtRest, G4bool wrappedIsAlongStep, G4bool wrappedIsPostStep,
82  G4String useThisName)
83  : G4VProcess( useThisName != "" ? useThisName : "biasWrapper("+wrappedProcess->GetProcessName()+")",
84  wrappedProcess->GetProcessType()),
85  fCurrentTrack ( nullptr ),
86  fPreviousStepSize (-1.0), fCurrentMinimumStep( -1.0 ), fProposedSafety ( -1.0),
87  fOccurenceBiasingOperation( nullptr ), fFinalStateBiasingOperation( nullptr ), fNonPhysicsBiasingOperation( nullptr ),
88  fPreviousOccurenceBiasingOperation( nullptr ), fPreviousFinalStateBiasingOperation( nullptr ), fPreviousNonPhysicsBiasingOperation( nullptr ),
89  fResetWrappedProcessInteractionLength ( false ),
90  fWrappedProcess ( wrappedProcess ),
91  fIsPhysicsBasedBiasing ( true ),
92  fWrappedProcessIsAtRest ( wrappedIsAtRest ),
93  fWrappedProcessIsAlong ( wrappedIsAlongStep ),
94  fWrappedProcessIsPost ( wrappedIsPostStep ),
95  fWrappedProcessPostStepGPIL ( -1.0 ),
96  fBiasingPostStepGPIL ( -1.0 ),
97  fWrappedProcessInteractionLength ( -1.0 ),
98  fWrappedProcessForceCondition ( NotForced ),
99  fBiasingForceCondition ( NotForced ),
100  fWrappedProcessAlongStepGPIL ( -1.0 ),
101  fBiasingAlongStepGPIL ( -1.0 ),
102  fWrappedProcessGPILSelection ( NotCandidateForSelection ),
103  fBiasingGPILSelection ( NotCandidateForSelection ),
104  fBiasingInteractionLaw ( nullptr ),
105  fPreviousBiasingInteractionLaw ( nullptr ),
106  fPhysicalInteractionLaw ( nullptr ),
107  fOccurenceBiasingParticleChange ( nullptr ),
108  fIamFirstGPIL ( false ),
109  fProcessManager ( nullptr ),
110  fSharedData ( nullptr )
111 {
112  for (G4int i = 0 ; i < 8 ; i++) fFirstLastFlags[i] = false;
113  fResetInteractionLaws.Put( true );
114  fCommonStart.Put(true);
115  fCommonEnd.Put(true);
116  fDoCommonConfigure.Put(true);
117 
119 
120  // -- create physical interaction law:
121  fPhysicalInteractionLaw = new G4InteractionLawPhysical("PhysicalInteractionLawFor("+GetProcessName()+")");
122  // -- instantiate particle change wrapper for occurrence biaising:
124  // -- instantiate a "do nothing" particle change:
126 }
127 
128 
129 
131 {
135 }
136 
137 
139 {
140  G4MapCache< const G4ProcessManager*,
142  if ( itr != G4BiasingProcessSharedData::fSharedDataMap.End( ) )
143  {
144  return (*itr).second;
145  }
146  else return 0;
147 }
148 
149 
151 {
162 
163  fPreviousStepSize = -1.0;
164 
166 
167  if ( fCommonStart.Get() )
168  {
169  fCommonStart.Put( false );// = false;
170  fCommonEnd .Put( true );// = true;
171 
172  fSharedData-> fCurrentBiasingOperator = 0;
174 
175  // -- §§ Add a "fSharedData->nStarting" here and outside bracket "fSharedData->nStarting++" and " if (fSharedData->nStarting) == fSharedData->(vector interface length)"
176  // -- §§ call to the loop "StartTracking" of operators"
177 
178  for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
180  }
181 }
182 
183 
185 {
189 
190  // -- Inform operators of end of tracking:
191  if ( fCommonEnd.Get() )
192  {
193  fCommonEnd .Put( false );// = false;
194  fCommonStart.Put( true );// = true;
195 
196  for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
198 
199  // -- §§ for above loop, do as in StartTracking.
200  }
201 }
202 
203 
204 
206  G4double previousStepSize,
208 {
209 
210  // ---------------------------------------------------------------------------------------------------
211  // -- The "biasing process master" takes care of updating the biasing operator, and for all biasing
212  // -- processes it invokes the PostStepGPIL of physical wrapped processes (anticipate stepping manager
213  // -- call ! ) to make all cross-sections updated with current step, and hence available before the
214  // -- first call to the biasing operator.
215  // ---------------------------------------------------------------------------------------------------
216  if ( fIamFirstGPIL )
217  {
218  // -- Update previous biasing operator, and assume the operator stays the same by
219  // -- default and that it is not left at the beginning of this step. These
220  // -- assumptions might be wrong if there is a volume change (in paralllel or
221  // -- mass geometries) in what case the flags will be updated.
223  fSharedData->fIsNewOperator = false;
225  // -- If new volume, either in mass or parallel geometries, get possible new biasing operator:
226  // -------------------------------------------------------------------------------------------
227  // -- Get biasing operator in parallel geometries:
228  G4bool firstStepInParallelVolume = false;
230  {
231  G4VBiasingOperator* newParallelOperator( nullptr );
232  G4bool firstStep = ( track.GetCurrentStepNumber() == 1 );
233  size_t iParallel = 0;
234  for ( auto wasLimiting : fSharedData->fParallelGeometriesLimiterProcess->GetWasLimiting() )
235  {
236  if ( firstStep || wasLimiting )
237  {
238  firstStepInParallelVolume = true;
239 
241  ->GetLogicalVolume() );
242  if ( newParallelOperator )
243  {
244  if ( tmpParallelOperator )
245  {
247  ed << " Several biasing operators are defined at the same place in parallel geometries ! Found:\n";
248  ed << " - `" << newParallelOperator->GetName() << "' and \n";
249  ed << " - `" << tmpParallelOperator->GetName() << "'.\n";
250  ed << " Keeping `" << newParallelOperator->GetName() << "'. Behavior not guaranteed ! Please consider having only one operator at a place. " << G4endl;
251  G4Exception(" G4BiasingProcessInterface::PostStepGetPhysicalInteractionLength(...)",
252  "BIAS.GEN.30",
253  JustWarning,
254  ed);
255  }
256  }
257  else newParallelOperator = tmpParallelOperator;
258  }
259  iParallel++;
260  }
261  fSharedData->fParallelGeometryOperator = newParallelOperator;
262  } // -- end of " if ( fSharedData->fParallelGeometriesLimiterProcess )"
263 
264  // -- Get biasing operator in mass geometry:
265  // -- [§§ Note : bug with this first step ? Does not work if previous step was concurrently limited with geometry. Might make use of safety at last point ?]
266  G4bool firstStepInVolume = ( (track.GetStep()->GetPreStepPoint()->GetStepStatus() == fGeomBoundary) || (track.GetCurrentStepNumber() == 1) );
267  // fSharedData->fIsNewOperator = false;
268  // fSharedData->fLeavingPreviousOperator = false;
269  if ( firstStepInVolume )
270  {
272  fSharedData->fMassGeometryOperator = newOperator;
273  if ( ( newOperator != nullptr ) && ( fSharedData->fParallelGeometryOperator != nullptr ) )
274  {
276  ed << " Biasing operators are defined at the same place in mass and parallel geometries ! Found:\n";
277  ed << " - `" << fSharedData->fParallelGeometryOperator->GetName() << "' in parallel geometry and \n";
278  ed << " - `" << newOperator->GetName() << "' in mass geometry.\n";
279  ed << " Keeping `" << fSharedData->fParallelGeometryOperator->GetName() << "'. Behavior not guaranteed ! Please consider having only one operator at a place. " << G4endl;
280  G4Exception(" G4BiasingProcessInterface::PostStepGetPhysicalInteractionLength(...)",
281  "BIAS.GEN.31",
282  JustWarning,
283  ed);
284  }
285  }
286 
287  // -- conclude the operator selection, giving priority to parallel geometry (as told in exception message BIAS.GEN.30):
288  if ( firstStepInVolume || firstStepInParallelVolume )
289  {
291  if ( newOperator == nullptr ) newOperator = fSharedData->fMassGeometryOperator;
292 
293  fSharedData->fCurrentBiasingOperator = newOperator ;
294 
295  if ( newOperator != fSharedData->fPreviousBiasingOperator )
296  {
298  fSharedData->fIsNewOperator = ( newOperator != nullptr );
299  }
300  }
301 
302 
303  // -- calls to wrapped process PostStepGPIL's:
304  // -------------------------------------------
305  // -- Each physics wrapper process has its
306  // -- fWrappedProcessPostStepGPIL ,
307  // -- fWrappedProcessForceCondition ,
308  // -- fWrappedProcessInteractionLength
309  // -- updated.
310  if ( fSharedData->fCurrentBiasingOperator != nullptr )
311  {
312  for ( size_t i = 0 ; i < (fSharedData->fPhysicsBiasingProcessInterfaces).size(); i++ )
313  (fSharedData->fPhysicsBiasingProcessInterfaces)[i]->InvokeWrappedProcessPostStepGPIL( track, previousStepSize, condition );
314  }
315  } // -- end of "if ( fIamFirstGPIL )"
316 
317 
318 
319  // -- Remember previous operator and proposed operations, if any, and reset:
320  // -------------------------------------------------------------------------
321  // -- remember only in case some biasing might be called
322  if ( ( fSharedData->fPreviousBiasingOperator != 0 ) ||
324  {
329  // -- reset:
334  // -- Physics PostStep and AlongStep GPIL
335  // fWrappedProcessPostStepGPIL : updated by InvokeWrappedProcessPostStepGPIL(...) above
337  // fWrappedProcessInteractionLength : updated by InvokeWrappedProcessPostStepGPIL(...) above; inverse of analog cross-section.
338  // fWrappedProcessForceCondition : updated by InvokeWrappedProcessPostStepGPIL(...) above
344  // -- for helper:
345  fPreviousStepSize = previousStepSize;
346  }
347 
348 
349  // -- previous step size value; it is switched to zero if resetting a wrapped process:
350  // -- (same trick used than in InvokedWrappedProcessPostStepGPIL )
351  G4double usedPreviousStepSize = previousStepSize;
352 
353  // ----------------------------------------------
354  // -- If leaving a biasing operator, let it know:
355  // ----------------------------------------------
357  {
358  (fSharedData->fPreviousBiasingOperator)->ExitingBiasing( &track, this );
359  // -- if no further biasing operator, reset process behavior to standard tracking:
361  {
364  {
365  // -- if the physics process has been under occurrence biasing, reset it:
367  {
370  // -- We set "previous step size" as 0.0, to let the process believe this is first step:
371  usedPreviousStepSize = 0.0;
372  }
373  }
374  }
375  }
376 
377 
378  // --------------------------------------------------------------
379  // -- no operator : analog tracking if physics-based, or nothing:
380  // --------------------------------------------------------------
382  {
383  // -- take note of the "usedPreviousStepSize" value:
384  if ( fIsPhysicsBasedBiasing ) return fWrappedProcess->PostStepGetPhysicalInteractionLength(track, usedPreviousStepSize, condition);
385  else
386  {
387  *condition = NotForced;
388  return DBL_MAX;
389  }
390  }
391 
392  // --------------------------------------------------
393  // -- A biasing operator exists. Proceed with
394  // -- treating non-physics and physics biasing cases:
395  //---------------------------------------------------
396 
397  // -- non-physics-based biasing case:
398  // ----------------------------------
399  if ( !fIsPhysicsBasedBiasing )
400  {
401  fNonPhysicsBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedNonPhysicsBiasingOperation( &track, this );
402  if ( fNonPhysicsBiasingOperation == 0 )
403  {
404  *condition = NotForced;
405  return DBL_MAX;
406  }
407  return fNonPhysicsBiasingOperation->DistanceToApplyOperation(&track, previousStepSize, condition);
408  }
409 
410 
411  // -- Physics based biasing case:
412  // ------------------------------
413  // -- Ask for possible GPIL biasing operation:
414  fOccurenceBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedOccurenceBiasingOperation( &track, this );
415 
416 
417  // -- no operation for occurrence biasing, analog GPIL returns the wrapped process GPIL and condition values
418  if ( fOccurenceBiasingOperation == 0 )
419  {
420  *condition = fWrappedProcessForceCondition;
422  }
423 
424  // -- A valid GPIL biasing operation has been proposed:
425  // -- 0) remember wrapped process will need to be reset on biasing exit, if particle survives:
427  // -- 1) update process interaction length for reference analog interaction law ( fWrappedProcessInteractionLength updated/collected above):
429  // -- 2) Collect biasing interaction law:
430  // -- The interaction law pointer is collected as a const pointer to the interaction law object.
431  // -- This interaction law will be kept under control of the biasing operation, which is the only
432  // -- entity that will change the state of the biasing interaction law.
433  // -- The force condition for biasing is asked at the same time, passing the analog one as default:
436  // -- 3) Ask operation to sample the biasing interaction law:
438 
439  // -- finish
440  *condition = fBiasingForceCondition;
441  return fBiasingPostStepGPIL;
442 
443 }
444 
445 
446 
448  const G4Step& step)
449 {
450  // ---------------------------------------
451  // -- case outside of volume with biasing:
452  // ---------------------------------------
453  if ( fSharedData->fCurrentBiasingOperator == 0 ) return fWrappedProcess->PostStepDoIt(track, step);
454 
455  // ----------------------------
456  // -- non-physics biasing case:
457  // ----------------------------
458  if ( !fIsPhysicsBasedBiasing )
459  {
460  G4VParticleChange* particleChange = fNonPhysicsBiasingOperation->GenerateBiasingFinalState( &track, &step );
461  (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC_NonPhysics, fNonPhysicsBiasingOperation, particleChange );
462  return particleChange;
463  }
464 
465  // -- physics biasing case:
466  // ------------------------
467  // -- It proceeds with the following logic:
468  // -- 1) Obtain the final state
469  // -- This final state may be analog or biased.
470  // -- The biased final state is obtained through a biasing operator
471  // -- returned by the operator.
472  // -- 2) The biased final state may be asked to be "force as it is"
473  // -- in what case the particle change is returned as is to the
474  // -- stepping.
475  // -- In all other cases (analog final state or biased final but
476  // -- not forced) the final state weight may be modified by the
477  // -- occurrence biasing, if such an occurrence biasing is at play.
478 
479  // -- Get final state, biased or analog:
480  G4VParticleChange* finalStateParticleChange;
482  fFinalStateBiasingOperation = (fSharedData->fCurrentBiasingOperator)->GetProposedFinalStateBiasingOperation( &track, this );
483  // -- Flag below is to force the biased generated particle change to be returned "as is" to the stepping, disregarding there
484  // -- was or not a occurrence biasing that would apply. Weight relevance under full responsibility of the biasing operation.
485  G4bool forceBiasedFinalState = false;
486  if ( fFinalStateBiasingOperation != 0 )
487  {
488  finalStateParticleChange = fFinalStateBiasingOperation->ApplyFinalStateBiasing( this, &track, &step, forceBiasedFinalState );
489  BAC = BAC_FinalState;
490  }
491  else
492  {
493  finalStateParticleChange = fWrappedProcess->PostStepDoIt(track, step);
494  BAC = BAC_None ;
495  }
496 
497  // -- if no occurrence biasing operation, we're done:
498  if ( fOccurenceBiasingOperation == 0 )
499  {
500  (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC, fFinalStateBiasingOperation, finalStateParticleChange );
501  return finalStateParticleChange;
502  }
503 
504  // -- if biased final state has been asked to be forced, we're done:
505  if ( forceBiasedFinalState )
506  {
507  (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC, fFinalStateBiasingOperation, finalStateParticleChange );
508  return finalStateParticleChange;
509  }
510 
511 
512  // -- If occurrence biasing, applies the occurrence biasing weight correction on top of final state (biased or not):
513  G4double weightForInteraction = 1.0;
514  if ( !fBiasingInteractionLaw->IsSingular() ) weightForInteraction =
516  fBiasingInteractionLaw ->ComputeEffectiveCrossSectionAt(step.GetStepLength());
517  else
518  {
519  // -- at this point effective XS can only be infinite, if not, there is a logic problem
521  {
523  ed << "Internal inconsistency in cross-section handling. Please report !" << G4endl;
524  G4Exception(" G4BiasingProcessInterface::PostStepDoIt(...)",
525  "BIAS.GEN.02",
526  JustWarning,
527  ed);
528  // -- if XS is infinite, weight is zero (and will stay zero), but we'll do differently.
529  // -- Should foresee in addition something to remember that in case of singular
530  // -- distribution, weight can only be partly calculated
531  }
532  }
533 
534  if ( weightForInteraction <= 0. )
535  {
537  ed << " Negative interaction weight : w_I = "
538  << weightForInteraction <<
541  " step length = " << step.GetStepLength() <<
542  " Interaction law = `" << fBiasingInteractionLaw << "'" <<
543  G4endl;
544  G4Exception(" G4BiasingProcessInterface::PostStepDoIt(...)",
545  "BIAS.GEN.03",
546  JustWarning,
547  ed);
548 
549  }
550 
551  (fSharedData->fCurrentBiasingOperator)->ReportOperationApplied( this, BAC,
552  fOccurenceBiasingOperation, weightForInteraction,
553  fFinalStateBiasingOperation, finalStateParticleChange );
554 
555 
558  fOccurenceBiasingParticleChange->SetWrappedParticleChange( finalStateParticleChange );
559  fOccurenceBiasingParticleChange->ProposeTrackStatus( finalStateParticleChange->GetTrackStatus() );
560  fOccurenceBiasingParticleChange->StealSecondaries(); // -- this also makes weightForInteraction applied to secondaries stolen
561 
562  // -- finish:
564 
565 }
566 
567 
568 // -- AlongStep methods:
570  G4double previousStepSize,
571  G4double currentMinimumStep,
572  G4double& proposedSafety,
573  G4GPILSelection* selection)
574 {
575 
576  // -- for helper methods:
577  fCurrentMinimumStep = currentMinimumStep;
578  fProposedSafety = proposedSafety;
579 
580 
581  // -- initialization default case:
583  *selection = NotCandidateForSelection;
584  // ---------------------------------------
585  // -- case outside of volume with biasing:
586  // ---------------------------------------
588  {
591  previousStepSize,
592  currentMinimumStep,
593  proposedSafety,
594  selection);
596  }
597 
598  // --------------------------------------------------------------------
599  // -- non-physics based biasing: no along operation expected (for now):
600  // --------------------------------------------------------------------
602 
603  // ----------------------
604  // -- physics-based case:
605  // ----------------------
606  if ( fOccurenceBiasingOperation == 0 )
607  {
610  previousStepSize,
611  currentMinimumStep,
612  proposedSafety,
613  selection);
615  }
616 
617 
618  // -----------------------------------------------------------
619  // -- From here we have an valid occurrence biasing operation:
620  // -----------------------------------------------------------
621  // -- Give operation opportunity to shorten step proposed by physics process:
623  G4double minimumStep = fBiasingAlongStepGPIL < currentMinimumStep ? fBiasingAlongStepGPIL : currentMinimumStep ;
624  // -- wrapped process is called with minimum step ( <= currentMinimumStep passed ) : an along process can not
625  // -- have its operation stretched over what it expects:
627  {
629  previousStepSize,
630  minimumStep,
631  proposedSafety,
632  selection);
633  fWrappedProcessGPILSelection = *selection;
635  }
636  else
637  {
640  }
641 
642  *selection = fBiasingGPILSelection;
643 
645 
646 }
647 
649  const G4Step& step)
650 {
651 
652  // ---------------------------------------
653  // -- case outside of volume with biasing:
654  // ---------------------------------------
656  {
657  if ( fWrappedProcessIsAlong ) return fWrappedProcess->AlongStepDoIt(track, step);
658  else
659  {
661  return fDummyParticleChange;
662  }
663  }
664 
665  // -----------------------------------
666  // -- case inside volume with biasing:
667  // -----------------------------------
669  else
670  {
673  }
674  G4double weightForNonInteraction (1.0);
675  if ( fBiasingInteractionLaw != 0 )
676  {
677  weightForNonInteraction =
679  fBiasingInteractionLaw ->ComputeNonInteractionProbabilityAt(step.GetStepLength());
680 
681  fOccurenceBiasingOperation->AlongMoveBy( this, &step, weightForNonInteraction );
682 
683  if ( weightForNonInteraction <= 0. )
684  {
686  ed << " Negative non interaction weight : w_NI = " << weightForNonInteraction <<
688  " p_NI(bias) = " << fBiasingInteractionLaw ->ComputeNonInteractionProbabilityAt(step.GetStepLength()) <<
689  " step length = " << step.GetStepLength() <<
690  " biasing interaction law = `" << fBiasingInteractionLaw->GetName() << "'" << G4endl;
691  G4Exception(" G4BiasingProcessInterface::AlongStepDoIt(...)",
692  "BIAS.GEN.04",
693  JustWarning,
694  ed);
695  }
696 
697  }
698 
700 
702 
703 }
704 
705 // -- AtRest methods
708 {
709  return fWrappedProcess->AtRestGetPhysicalInteractionLength(track, condition);
710 }
712  const G4Step& step)
713 {
714  return fWrappedProcess->AtRestDoIt(track, step);
715 }
716 
717 
719 {
720  if ( fWrappedProcess != 0 ) return fWrappedProcess->IsApplicable(pd);
721  else return true;
722 }
723 
724 
726 {
727  // -- Master for this process:
729  // -- Master for wrapped process:
730  if ( fWrappedProcess != 0 )
731  {
732  const G4BiasingProcessInterface* thisWrapperMaster = (const G4BiasingProcessInterface *)GetMasterProcess();
733  // -- paranoia check: (?)
734  G4VProcess* wrappedMaster = 0;
735  wrappedMaster = thisWrapperMaster->GetWrappedProcess();
736  fWrappedProcess->SetMasterProcess( wrappedMaster );
737  }
738 }
739 
740 
742 {
743  // -- Sequential mode : called second (after PreparePhysicsTable(..))
744  // -- MT mode : called second (after PreparePhysicsTable(..)) by master thread.
745  // -- Corresponding process instance not used then by tracking.
746  // -- PreparePhysicsTable(...) has been called first for all processes,
747  // -- so the first/last flags and G4BiasingProcessInterface vector of processes have
748  // -- been properly setup, fIamFirstGPIL is valid.
749  if ( fWrappedProcess != 0 )
750  {
752  }
753 
754  if ( fIamFirstGPIL )
755  {
756  // -- Re-order vector of processes to match that of the GPIL
757  // -- (made for fIamFirstGPIL, but important is to have it made once):
759  // -- Let operators to configure themselves for the master thread or for sequential mode.
760  // -- Intended here is in particular the registration to physics model catalog.
761  // -- The fDoCommonConfigure is to ensure that this Configure is made by only one process (othewise each first process makes the call):
762  if ( fDoCommonConfigure.Get() )
763  {
764  for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
766  fDoCommonConfigure.Put(false);
767  }
768 
769  }
770 }
771 
772 
774 {
775  // -- Sequential mode : called first (before BuildPhysicsTable(..))
776  // -- MT mode : called first (before BuildPhysicsTable(..)) by master thread.
777  // -- Corresponding process instance not used then by tracking.
778  // -- Let process finding its first/last position in the process manager:
780  if ( fWrappedProcess != 0 )
781  {
783  }
784 }
785 
786 
788 {
789  if ( fWrappedProcess != 0 ) return fWrappedProcess->StorePhysicsTable(pd, s, f);
790  else return false;
791 }
792 
793 
795 {
796  if ( fWrappedProcess != 0 ) return fWrappedProcess->RetrievePhysicsTable(pd, s, f);
797  else return false;
798 }
799 
800 
802 {
805 
806  // -- initialize fSharedData pointer:
808  {
811  }
813  // -- augment list of co-operating processes:
814  fSharedData-> fBiasingProcessInterfaces.push_back( this );
815  fSharedData-> fPublicBiasingProcessInterfaces.push_back( this );
816  if ( fIsPhysicsBasedBiasing )
817  {
818  fSharedData-> fPhysicsBiasingProcessInterfaces.push_back( this );
819  fSharedData-> fPublicPhysicsBiasingProcessInterfaces.push_back( this );
820  }
821  else
822  {
823  fSharedData-> fNonPhysicsBiasingProcessInterfaces.push_back( this );
824  fSharedData-> fPublicNonPhysicsBiasingProcessInterfaces.push_back( this );
825  }
826  // -- remember process manager:
827  fProcessManager = mgr;
828 }
829 
830 
832 {
833  if ( fWrappedProcess != 0 ) return fWrappedProcess->GetProcessManager();
834  else return G4VProcess::GetProcessManager();
835 }
836 
837 
839 {
840  // -- Sequential mode : not called
841  // -- MT mode : called after PrepareWorkerPhysicsTable(..)
842  // -- PrepareWorkerPhysicsTable(...) has been called first for all processes,
843  // -- so the first/last flags and G4BiasingProcessInterface vector of processes have
844  // -- been properly setup, fIamFirstGPIL is valid.
845  if ( fWrappedProcess != 0 )
846  {
848  }
849 
850  if ( fIamFirstGPIL )
851  {
852  // -- Re-order vector of processes to match that of the GPIL
853  // -- (made for fIamFirstGPIL, but important is to have it made once):
855  // -- Let operators to configure themselves for the worker thread, if needed.
856  // -- Registration to physics model catalog **IS NOT** to be made here, but in Configure().
857  // -- The fDoCommonConfigure is to ensure that this Configure is made by only one process (othewise each first process makes the call):
858  if ( fDoCommonConfigure.Get() )
859  {
860  for ( size_t optr = 0 ; optr < ( G4VBiasingOperator::GetBiasingOperators() ).size() ; optr ++)
861  ( G4VBiasingOperator::GetBiasingOperators() )[optr]->ConfigureForWorker( );
862  fDoCommonConfigure.Put(false);
863  }
864  }
865 }
866 
867 
869 {
870  // -- Sequential mode : not called
871  // -- MT mode : called first, before BuildWorkerPhysicsTable(..)
872  // -- Let process finding its first/last position in the process manager:
874 
875  if ( fWrappedProcess != 0 )
876  {
878  }
879 }
880 
881 
883 {
885 }
886 
887 
889 {
890  G4int iPhys = ( physOnly ) ? 1 : 0;
891  return fFirstLastFlags[IdxFirstLast( 1, 1, iPhys)];
892 }
893 
894 
896 {
897  G4int iPhys = ( physOnly ) ? 1 : 0;
898  return fFirstLastFlags[IdxFirstLast( 0, 1, iPhys)];
899 }
900 
901 
903 {
904  G4int iPhys = ( physOnly ) ? 1 : 0;
905  return fFirstLastFlags[IdxFirstLast( 1, 0, iPhys)];
906 }
907 
908 
910 {
911  G4int iPhys = ( physOnly ) ? 1 : 0;
912  return fFirstLastFlags[IdxFirstLast( 0, 0, iPhys)];
913 }
914 
915 
917 {
918  G4bool isFirst = true;
920  G4int thisIdx(-1);
921  for (std::size_t i = 0; i < pv->size(); ++i ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
922  if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
923  for ( std::size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); ++i )
924  {
925  if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
926  {
927  G4int thatIdx(-1);
928  for (std::size_t j = 0; j < pv->size(); ++j ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
929  if ( thatIdx >= 0 ) // -- to ignore pure along processes
930  {
931  if ( thisIdx > thatIdx )
932  {
933  isFirst = false;
934  break;
935  }
936  }
937  }
938  }
939  return isFirst;
940 }
941 
942 
944 {
945  G4bool isLast = true;
947  G4int thisIdx(-1);
948  for (std::size_t i = 0; i < pv->size(); ++i ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
949  if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
950  for ( std::size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); ++i )
951  {
952  if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
953  {
954  G4int thatIdx(-1);
955  for (std::size_t j = 0; j < pv->size(); ++j ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
956  if ( thatIdx >= 0 ) // -- to ignore pure along processes
957  {
958  if ( thisIdx < thatIdx )
959  {
960  isLast = false;
961  break;
962  }
963  }
964  }
965  }
966  return isLast;
967 }
968 
969 
971 {
972  G4bool isFirst = true;
974  G4int thisIdx(-1);
975  for (std::size_t i = 0; i < pv->size(); ++i ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
976  if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
977  for ( std::size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); ++i )
978  {
979  if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
980  {
981  G4int thatIdx(-1);
982  for (std::size_t j = 0; j < pv->size(); ++j ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
983  if ( thatIdx >= 0 ) // -- to ignore pure along processes
984  {
985  if ( thisIdx > thatIdx )
986  {
987  isFirst = false;
988  break;
989  }
990  }
991  }
992  }
993  return isFirst;
994 }
995 
996 
998 {
999  G4bool isLast = true;
1001  G4int thisIdx(-1);
1002  for (std::size_t i = 0; i < pv->size(); ++i ) if ( (*pv)(i) == this ) { thisIdx = i; break; }
1003  if ( thisIdx < 0 ) return false; // -- to ignore pure along processes
1004  for ( std::size_t i = 0; i < (fSharedData->fBiasingProcessInterfaces).size(); ++i )
1005  {
1006  if ( (fSharedData->fBiasingProcessInterfaces)[i]->fIsPhysicsBasedBiasing || !physOnly )
1007  {
1008  G4int thatIdx(-1);
1009  for (std::size_t j = 0; j < pv->size(); ++j ) if ( (*pv)(j) == (fSharedData->fBiasingProcessInterfaces)[i] ) { thatIdx = j; break; }
1010  if ( thatIdx >= 0 ) // -- to ignore pure along processes
1011  {
1012  if ( thisIdx < thatIdx )
1013  {
1014  isLast = false;
1015  break;
1016  }
1017  }
1018  }
1019  }
1020  return isLast;
1021 }
1022 
1023 
1025 {
1026  for ( G4int iPhys = 0; iPhys < 2; iPhys++ )
1027  {
1028  G4bool physOnly = ( iPhys == 1 );
1029  fFirstLastFlags[IdxFirstLast( 1, 1, iPhys)] = IsFirstPostStepGPILInterface(physOnly);
1030  fFirstLastFlags[IdxFirstLast( 0, 1, iPhys)] = IsLastPostStepGPILInterface(physOnly);
1031  fFirstLastFlags[IdxFirstLast( 1, 0, iPhys)] = IsFirstPostStepDoItInterface(physOnly);
1032  fFirstLastFlags[IdxFirstLast( 0, 0, iPhys)] = IsLastPostStepDoItInterface(physOnly);
1033  }
1034 
1035  // -- for itself, for optimization:
1037 }
1038 
1039 
1041 {
1046 }
1047 
1048 
1050  G4double previousStepSize,
1052 {
1053  G4double usedPreviousStepSize = previousStepSize;
1054  // -- if the physics process has been under occurrence biasing in the previous step
1055  // -- we reset it, as we don't know if it will be biased again or not in this
1056  // -- step. The pity is that PostStepGPIL and interaction length (cross-section)
1057  // -- calculations are done both in the PostStepGPIL of the process, while here we
1058  // -- are just interested in the calculation of the cross-section. This is a pity
1059  // -- as this forces to re-generated a random number for nothing.
1061  {
1064  // -- We set "previous step size" as 0.0, to let the process believe this is first step:
1065  usedPreviousStepSize = 0.0;
1066  }
1067  // -- GPIL response:
1068  fWrappedProcessPostStepGPIL = fWrappedProcess->PostStepGetPhysicalInteractionLength(track, usedPreviousStepSize, condition);
1070  // -- and (inverse) cross-section:
1072 }
1073 
1074 
1076 {
1077  // -- re-order vector of processes to match that of the GPIL:
1078  std::vector < G4BiasingProcessInterface* > tmpProcess ( fSharedData->fBiasingProcessInterfaces );
1079  ( fSharedData -> fBiasingProcessInterfaces ) . clear();
1080  ( fSharedData -> fPhysicsBiasingProcessInterfaces ) . clear();
1081  ( fSharedData -> fNonPhysicsBiasingProcessInterfaces ) . clear();
1082  ( fSharedData -> fPublicBiasingProcessInterfaces ) . clear();
1083  ( fSharedData -> fPublicPhysicsBiasingProcessInterfaces ) . clear();
1084  ( fSharedData -> fPublicNonPhysicsBiasingProcessInterfaces ) . clear();
1085 
1087  for (std::size_t i = 0; i < pv->size(); ++i )
1088  {
1089  for ( std::size_t j = 0; j < tmpProcess.size(); ++j )
1090  {
1091  if ( (*pv)(i) == tmpProcess[j] )
1092  {
1093  ( fSharedData -> fBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1094  ( fSharedData -> fPublicBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1095  if ( tmpProcess[j] -> fIsPhysicsBasedBiasing )
1096  {
1097  ( fSharedData -> fPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1098  ( fSharedData -> fPublicPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1099  }
1100  else
1101  {
1102  ( fSharedData -> fNonPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1103  ( fSharedData -> fPublicNonPhysicsBiasingProcessInterfaces ) . push_back( tmpProcess[j] );
1104  }
1105  break;
1106  }
1107  }
1108  }
1109 }