Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/biasing/generic/src/G4BiasingProcessInterface.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/biasing/generic/src/G4BiasingProcessInterface.cc (Version 11.3.0) and /processes/biasing/generic/src/G4BiasingProcessInterface.cc (Version 9.3.p1)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 // G4BiasingProcessInterface                      
 27 // -------------------------------------------    
 28                                                   
 29 #include "G4BiasingProcessInterface.hh"           
 30 #include "G4VBiasingOperator.hh"                  
 31 #include "G4VBiasingOperation.hh"                 
 32 #include "G4ParticleChangeForOccurenceBiasing.    
 33 #include "G4ParticleChangeForNothing.hh"          
 34 #include "G4VBiasingInteractionLaw.hh"            
 35 #include "G4InteractionLawPhysical.hh"            
 36 #include "G4ProcessManager.hh"                    
 37 #include "G4BiasingAppliedCase.hh"                
 38 #include "G4ParallelGeometriesLimiterProcess.h    
 39                                                   
 40 G4Cache<G4bool> G4BiasingProcessInterface::fRe    
 41 G4Cache<G4bool> G4BiasingProcessInterface::fCo    
 42 G4Cache<G4bool> G4BiasingProcessInterface::fCo    
 43 G4Cache<G4bool> G4BiasingProcessInterface::fDo    
 44                                                   
 45 G4BiasingProcessInterface::G4BiasingProcessInt    
 46   : G4VProcess( name ),                           
 47     fResetWrappedProcessInteractionLength( tru    
 48 {                                                 
 49   for (G4int i = 0 ; i < 8 ; i++)  fFirstLastF    
 50   fResetInteractionLaws.Put( true );              
 51   fCommonStart         .Put( true );              
 52   fCommonEnd           .Put( true );              
 53   fDoCommonConfigure   .Put( true );              
 54 }                                                 
 55                                                   
 56                                                   
 57 G4BiasingProcessInterface::                       
 58 G4BiasingProcessInterface(G4VProcess* wrappedP    
 59                           G4bool wrappedIsAtRe    
 60                           G4bool wrappedIsAlon    
 61                           G4bool wrappedIsPost    
 62                           G4String useThisName    
 63   : G4VProcess(useThisName != ""                  
 64                  ? useThisName                    
 65                  : "biasWrapper("+wrappedProce    
 66                wrappedProcess->GetProcessType(    
 67     fWrappedProcess( wrappedProcess ),            
 68     fIsPhysicsBasedBiasing( true ),               
 69     fWrappedProcessIsAtRest( wrappedIsAtRest )    
 70     fWrappedProcessIsAlong( wrappedIsAlongStep    
 71     fWrappedProcessIsPost( wrappedIsPostStep )    
 72 {                                                 
 73   for (G4int i = 0 ; i < 8 ; ++i)                 
 74     fFirstLastFlags[i] = false;                   
 75   fResetInteractionLaws.Put( true );              
 76   fCommonStart.Put(true);                         
 77   fCommonEnd.Put(true);                           
 78   fDoCommonConfigure.Put(true);                   
 79                                                   
 80   SetProcessSubType(fWrappedProcess->GetProces    
 81                                                   
 82   // -- create physical interaction law:          
 83   fPhysicalInteractionLaw = new G4InteractionL    
 84   // -- instantiate particle change wrapper fo    
 85   fOccurenceBiasingParticleChange = new G4Part    
 86   // -- instantiate a "do nothing" particle ch    
 87   fDummyParticleChange = new G4ParticleChangeF    
 88 }                                                 
 89                                                   
 90 G4BiasingProcessInterface::~G4BiasingProcessIn    
 91 {                                                 
 92   delete fPhysicalInteractionLaw;                 
 93   delete fOccurenceBiasingParticleChange;         
 94   delete fDummyParticleChange;                    
 95 }                                                 
 96                                                   
 97 const G4BiasingProcessSharedData*                 
 98 G4BiasingProcessInterface::GetSharedData( cons    
 99 {                                                 
100   const auto & itr = G4BiasingProcessSharedDat    
101   if ( itr !=  G4BiasingProcessSharedData::fSh    
102   {                                               
103     return (*itr).second;                         
104   }                                               
105   else return nullptr;                            
106 }                                                 
107                                                   
108                                                   
109 void G4BiasingProcessInterface::StartTracking(    
110 {                                                 
111   fCurrentTrack = track;                          
112   if ( fIsPhysicsBasedBiasing ) fWrappedProces    
113   fOccurenceBiasingOperation          = nullpt    
114   fPreviousOccurenceBiasingOperation  = nullpt    
115   fFinalStateBiasingOperation         = nullpt    
116   fPreviousFinalStateBiasingOperation = nullpt    
117   fNonPhysicsBiasingOperation         = nullpt    
118   fPreviousNonPhysicsBiasingOperation = nullpt    
119   fBiasingInteractionLaw              = nullpt    
120   fPreviousBiasingInteractionLaw      = nullpt    
121                                                   
122   fPreviousStepSize = -1.0;                       
123                                                   
124   fResetWrappedProcessInteractionLength = fals    
125                                                   
126   if ( fCommonStart.Get() )                       
127   {                                               
128     fCommonStart.Put( false );// = false;         
129     fCommonEnd.Put( true  );// = true;            
130                                                   
131     fSharedData->fCurrentBiasingOperator = nul    
132     fSharedData->fPreviousBiasingOperator = nu    
133                                                   
134     // -- Add a  "fSharedData->nStarting" here    
135     // -- call to the loop "StartTracking" of     
136                                                   
137     for (std::size_t optr=0 ; optr<(G4VBiasing    
138     {                                             
139       (G4VBiasingOperator::GetBiasingOperators    
140     }                                             
141   }                                               
142 }                                                 
143                                                   
144 void G4BiasingProcessInterface::EndTracking()     
145 {                                                 
146   if ( fIsPhysicsBasedBiasing )                   
147     fWrappedProcess->EndTracking();               
148   if ( fSharedData->fCurrentBiasingOperator)      
149     (fSharedData->fCurrentBiasingOperator)->Ex    
150   fBiasingInteractionLaw = nullptr;               
151                                                   
152   // -- Inform operators of end of tracking:      
153   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     {                                             
160       (G4VBiasingOperator::GetBiasingOperators    
161     }                                             
162     // -- for above loop, do as in StartTracki    
163   }                                               
164 }                                                 
165                                                   
166 G4double G4BiasingProcessInterface::              
167 PostStepGetPhysicalInteractionLength( const G4    
168                                       G4double    
169                                       G4ForceC    
170 {                                                 
171   // -----------------------------------------    
172   // -- The "biasing process master" takes car    
173   // -- processes it invokes the PostStepGPIL     
174   // -- call ! ) to make all cross-sections up    
175   // -- first call to the biasing operator.       
176   // -----------------------------------------    
177   if ( fIamFirstGPIL )                            
178   {                                               
179     // -- Update previous biasing operator, an    
180     // -- default and that it is not left at t    
181     // -- assumptions might be wrong if there     
182     // -- mass geometries) in what case the fl    
183     fSharedData->fPreviousBiasingOperator = fS    
184     fSharedData->fIsNewOperator           = fa    
185     fSharedData->fLeavingPreviousOperator = fa    
186     // -- If new volume, either in mass or par    
187     // ---------------------------------------    
188     // -- Get biasing operator in parallel geo    
189     G4bool firstStepInParallelVolume = false;     
190     if ( fSharedData->fParallelGeometriesLimit    
191     {                                             
192       G4VBiasingOperator* newParallelOperator(    
193       G4bool firstStep = ( track.GetCurrentSte    
194       std::size_t iParallel = 0;                  
195       for ( auto wasLimiting : fSharedData->fP    
196       {                                           
197         if ( firstStep || wasLimiting )           
198         {                                         
199           firstStepInParallelVolume = true;       
200           auto tmpParallelOperator = G4VBiasin    
201              GetBiasingOperator((fSharedData->    
202              ->GetCurrentVolumes()[iParallel])    
203              ->GetLogicalVolume());               
204           if ( newParallelOperator )              
205           {                                       
206             if ( tmpParallelOperator )            
207             {                                     
208               G4ExceptionDescription ed;          
209               ed << " Several biasing operator    
210                  << " in parallel geometries !    
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     }                                             
280   } // -- end of "if ( fIamFirstGPIL )"           
281                                                   
282   // -- Remember previous operator and propose    
283   // -----------------------------------------    
284   // -- remember only in case some biasing mig    
285   if ( ( fSharedData->fPreviousBiasingOperator    
286        ( fSharedData->fCurrentBiasingOperator     
287   {                                               
288     fPreviousOccurenceBiasingOperation  = fOcc    
289     fPreviousFinalStateBiasingOperation = fFin    
290     fPreviousNonPhysicsBiasingOperation = fNon    
291     fPreviousBiasingInteractionLaw      = fBia    
292     // -- reset:                                  
293     fOccurenceBiasingOperation          = null    
294     fFinalStateBiasingOperation         = null    
295     fNonPhysicsBiasingOperation         = null    
296     fBiasingInteractionLaw              = null    
297     // -- Physics PostStep and AlongStep GPIL     
298     // fWrappedProcessPostStepGPIL      : upda    
299     fBiasingPostStepGPIL                = DBL_    
300     // fWrappedProcessInteractionLength : upda    
301     // fWrappedProcessForceCondition    : upda    
302     fBiasingForceCondition              = NotF    
303     fWrappedProcessAlongStepGPIL        = DBL_    
304     fBiasingAlongStepGPIL               = DBL_    
305     fWrappedProcessGPILSelection        = NotC    
306     fBiasingGPILSelection               = NotC    
307     // -- for helper:                             
308     fPreviousStepSize                   = prev    
309   }                                               
310                                                   
311   // -- previous step size value; it is switch    
312   // -- (same trick used than in InvokedWrappe    
313   G4double usedPreviousStepSize = previousStep    
314                                                   
315   // -----------------------------------------    
316   // -- If leaving a biasing operator, let it     
317   // -----------------------------------------    
318   if ( fSharedData->fLeavingPreviousOperator )    
319   {                                               
320     (fSharedData->fPreviousBiasingOperator)->E    
321     // -- if no further biasing operator, rese    
322     if ( fSharedData->fCurrentBiasingOperator     
323     {                                             
324       ResetForUnbiasedTracking();                 
325       if ( fIsPhysicsBasedBiasing )               
326       {                                           
327         // -- if the physics process has been     
328         if ( fResetWrappedProcessInteractionLe    
329         {                                         
330           fResetWrappedProcessInteractionLengt    
331           fWrappedProcess->ResetNumberOfIntera    
332           // -- We set "previous step size" as    
333           usedPreviousStepSize = 0.0;             
334         }                                         
335       }                                           
336     }                                             
337   }                                               
338                                                   
339   // -----------------------------------------    
340   // -- no operator : analog tracking if physi    
341   // -----------------------------------------    
342   if ( fSharedData->fCurrentBiasingOperator ==    
343   {                                               
344     // -- take note of the "usedPreviousStepSi    
345     if ( fIsPhysicsBasedBiasing )                 
346     {                                             
347       return fWrappedProcess->PostStepGetPhysi    
348     }                                             
349     else                                          
350     {                                             
351       *condition = NotForced;                     
352       return DBL_MAX;                             
353     }                                             
354   }                                               
355                                                   
356   // -----------------------------------------    
357   // -- A biasing operator exists. Proceed wit    
358   // -- treating non-physics and physics biasi    
359   //------------------------------------------    
360                                                   
361   // -- non-physics-based biasing case:           
362   // ----------------------------------           
363   if ( !fIsPhysicsBasedBiasing )                  
364   {                                               
365     fNonPhysicsBiasingOperation = (fSharedData    
366     if ( fNonPhysicsBiasingOperation == nullpt    
367     {                                             
368       *condition = NotForced;                     
369       return DBL_MAX;                             
370     }                                             
371     return fNonPhysicsBiasingOperation->Distan    
372   }                                               
373                                                   
374   // -- Physics based biasing case:               
375   // ------------------------------               
376   // -- Ask for possible GPIL biasing operatio    
377   fOccurenceBiasingOperation = (fSharedData->f    
378                                                   
379   // -- no operation for occurrence biasing, a    
380   if ( fOccurenceBiasingOperation == nullptr )    
381   {                                               
382     *condition = fWrappedProcessForceCondition    
383     return fWrappedProcessPostStepGPIL;           
384   }                                               
385                                                   
386   // -- A valid GPIL biasing operation has bee    
387   // -- 0) remember wrapped process will need     
388   fResetWrappedProcessInteractionLength = true    
389   // -- 1) update process interaction length f    
390   fPhysicalInteractionLaw->SetPhysicalCrossSec    
391   // -- 2) Collect biasing interaction law:       
392   // --    The interaction law pointer is coll    
393   // --    This interaction law will be kept u    
394   // --    entity that will change the state o    
395   // --    The force condition for biasing is     
396   fBiasingForceCondition = fWrappedProcessForc    
397   fBiasingInteractionLaw = fOccurenceBiasingOp    
398   // -- 3) Ask operation to sample the biasing    
399   fBiasingPostStepGPIL = fBiasingInteractionLa    
400                                                   
401   // -- finish                                    
402   *condition = fBiasingForceCondition;            
403   return fBiasingPostStepGPIL;                    
404 }                                                 
405                                                   
406 G4VParticleChange* G4BiasingProcessInterface::    
407                                                   
408 {                                                 
409   // ---------------------------------------      
410   // -- case outside of volume with biasing:      
411   // ---------------------------------------      
412   if ( fSharedData->fCurrentBiasingOperator ==    
413     return fWrappedProcess->PostStepDoIt(track    
414                                                   
415   // ----------------------------                 
416   // -- non-physics biasing case:                 
417   // ----------------------------                 
418   if ( !fIsPhysicsBasedBiasing )                  
419   {                                               
420     G4VParticleChange* particleChange = fNonPh    
421     (fSharedData->fCurrentBiasingOperator)->Re    
422     return particleChange;                        
423   }                                               
424                                                   
425   // -- physics biasing case:                     
426   // ------------------------                     
427   // -- It proceeds with the following logic:     
428   // -- 1) Obtain the final state                 
429   // --    This final state may be analog or b    
430   // --    The biased final state is obtained     
431   // --    returned by the operator.              
432   // -- 2) The biased final state may be asked    
433   // --    in what case the particle change is    
434   // --    stepping.                              
435   // --    In all other cases (analog final st    
436   // --    not forced) the final state weight     
437   // --    occurrence biasing, if such an occu    
438                                                   
439   // -- Get final state, biased or analog:        
440   G4VParticleChange* finalStateParticleChange;    
441   G4BiasingAppliedCase BAC;                       
442   fFinalStateBiasingOperation = (fSharedData->    
443   // -- Flag below is to force the biased gene    
444   // -- was or not a occurrence biasing that w    
445   G4bool forceBiasedFinalState = false;           
446   if ( fFinalStateBiasingOperation != nullptr     
447   {                                               
448     finalStateParticleChange = fFinalStateBias    
449     BAC = BAC_FinalState;                         
450   }                                               
451   else                                            
452   {                                               
453     finalStateParticleChange = fWrappedProcess    
454     BAC =  BAC_None ;                             
455   }                                               
456                                                   
457   // -- if no occurrence biasing operation, we    
458   if ( fOccurenceBiasingOperation == nullptr )    
459   {                                               
460     (fSharedData->fCurrentBiasingOperator)->Re    
461     return finalStateParticleChange;              
462   }                                               
463                                                   
464   // -- if biased final state has been asked t    
465   if ( forceBiasedFinalState )                    
466   {                                               
467     (fSharedData->fCurrentBiasingOperator)->Re    
468     return finalStateParticleChange;              
469   }                                               
470                                                   
471   // -- If occurrence biasing, applies the occ    
472   G4double weightForInteraction = 1.0;            
473   if ( !fBiasingInteractionLaw->IsSingular() )    
474   {                                               
475     weightForInteraction = fPhysicalInteractio    
476                          / fBiasingInteraction    
477   }                                               
478   else                                            
479   {                                               
480     // -- at this point effective XS can only     
481     if ( !fBiasingInteractionLaw->IsEffectiveC    
482     {                                             
483       G4ExceptionDescription ed;                  
484       ed << "Internal inconsistency in cross-s    
485       G4Exception(" G4BiasingProcessInterface:    
486                   "BIAS.GEN.02", JustWarning,     
487       // -- if XS is infinite, weight is zero     
488       // -- Should foresee in addition somethi    
489       // -- distribution, weight can only be p    
490     }                                             
491   }                                               
492                                                   
493   if ( weightForInteraction <= 0. )               
494   {                                               
495     G4ExceptionDescription ed;                    
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                                                   
514   fOccurenceBiasingParticleChange->SetOccurenc    
515   fOccurenceBiasingParticleChange->SetSecondar    
516   fOccurenceBiasingParticleChange->SetWrappedP    
517   fOccurenceBiasingParticleChange->ProposeTrac    
518   fOccurenceBiasingParticleChange->StealSecond    
519                                                   
520   // -- finish:                                   
521   return fOccurenceBiasingParticleChange;         
522 }                                                 
523                                                   
524 // -- AlongStep methods:                          
525 G4double G4BiasingProcessInterface::              
526 AlongStepGetPhysicalInteractionLength(const G4    
527                                             G4    
528                                             G4    
529                                             G4    
530                                             G4    
531 {                                                 
532   // -- for helper methods:                       
533   fCurrentMinimumStep = currentMinimumStep;       
534   fProposedSafety = proposedSafety;               
535                                                   
536   // -- initialization default case:              
537   fWrappedProcessAlongStepGPIL = DBL_MAX;         
538   *selection                   = NotCandidateF    
539   // ---------------------------------------      
540   // -- case outside of volume with biasing:      
541   // ---------------------------------------      
542   if ( fSharedData->fCurrentBiasingOperator ==    
543   {                                               
544     if ( fWrappedProcessIsAlong )                 
545       fWrappedProcessAlongStepGPIL = fWrappedP    
546         ->AlongStepGetPhysicalInteractionLengt    
547                                                   
548                                                   
549     return fWrappedProcessAlongStepGPIL;          
550   }                                               
551                                                   
552   // -----------------------------------------    
553   // -- non-physics based biasing: no along op    
554   // -----------------------------------------    
555   if ( !fIsPhysicsBasedBiasing ) return fWrapp    
556                                                   
557   // ----------------------                       
558   // -- physics-based case:                       
559   // ----------------------                       
560   if ( fOccurenceBiasingOperation == nullptr )    
561   {                                               
562     if ( fWrappedProcessIsAlong )                 
563       fWrappedProcessAlongStepGPIL = fWrappedP    
564         ->AlongStepGetPhysicalInteractionLengt    
565                                                   
566                                                   
567     return fWrappedProcessAlongStepGPIL;          
568   }                                               
569                                                   
570   // -----------------------------------------    
571   // -- From here we have an valid occurrence     
572   // -----------------------------------------    
573   // -- Give operation opportunity to shorten     
574   fBiasingAlongStepGPIL = fOccurenceBiasingOpe    
575   G4double minimumStep  = fBiasingAlongStepGPI    
576                         ? fBiasingAlongStepGPI    
577   // -- wrapped process is called with minimum    
578   // -- have its operation stretched over what    
579   if ( fWrappedProcessIsAlong )                   
580   {                                               
581     fWrappedProcessAlongStepGPIL = fWrappedPro    
582       ->AlongStepGetPhysicalInteractionLength(    
583                                                   
584                                                   
585     fWrappedProcessGPILSelection = *selection;    
586     fBiasingGPILSelection = fOccurenceBiasingO    
587       ->ProposeGPILSelection( fWrappedProcessG    
588   }                                               
589   else                                            
590   {                                               
591     fBiasingGPILSelection = fOccurenceBiasingO    
592       ->ProposeGPILSelection( NotCandidateForS    
593     fWrappedProcessAlongStepGPIL = fBiasingAlo    
594   }                                               
595                                                   
596   *selection = fBiasingGPILSelection;             
597                                                   
598   return fWrappedProcessAlongStepGPIL;            
599 }                                                 
600                                                   
601 G4VParticleChange*                                
602 G4BiasingProcessInterface::AlongStepDoIt(const    
603                                          const    
604 {                                                 
605   // ---------------------------------------      
606   // -- case outside of volume with biasing:      
607   // ---------------------------------------      
608   if ( fSharedData->fCurrentBiasingOperator ==    
609   {                                               
610     if ( fWrappedProcessIsAlong )                 
611     {                                             
612       return fWrappedProcess->AlongStepDoIt(tr    
613     }                                             
614     else                                          
615     {                                             
616       fDummyParticleChange->Initialize( track     
617       return fDummyParticleChange;                
618     }                                             
619   }                                               
620                                                   
621   // -----------------------------------          
622   // -- case inside volume with biasing:          
623   // -----------------------------------          
624   if ( fWrappedProcessIsAlong )                   
625   {                                               
626     fOccurenceBiasingParticleChange               
627       ->SetWrappedParticleChange(fWrappedProce    
628   }                                               
629   else                                            
630   {                                               
631     fOccurenceBiasingParticleChange->SetWrappe    
632     fOccurenceBiasingParticleChange->ProposeTr    
633   }                                               
634   G4double weightForNonInteraction (1.0);         
635   if ( fBiasingInteractionLaw != nullptr )        
636   {                                               
637     weightForNonInteraction =                     
638       fPhysicalInteractionLaw->ComputeNonInter    
639       fBiasingInteractionLaw ->ComputeNonInter    
640                                                   
641     fOccurenceBiasingOperation->AlongMoveBy( t    
642                                                   
643     if ( weightForNonInteraction <= 0. )          
644     {                                             
645       G4ExceptionDescription ed;                  
646       ed << " Negative non interaction weight     
647             " p_NI(phys) = " <<  fPhysicalInte    
648             " p_NI(bias) = " <<  fBiasingInter    
649             " step length = "  <<  step.GetSte    
650             " biasing interaction law = `" <<     
651       G4Exception(" G4BiasingProcessInterface:    
652                   "BIAS.GEN.04", JustWarning,     
653     }                                             
654   }                                               
655                                                   
656   fOccurenceBiasingParticleChange                 
657     ->SetOccurenceWeightForNonInteraction( wei    
658                                                   
659   return fOccurenceBiasingParticleChange;         
660 }                                                 
661                                                   
662 // -- AtRest methods                              
663 G4double G4BiasingProcessInterface::              
664 AtRestGetPhysicalInteractionLength(const G4Tra    
665                                    G4ForceCond    
666 {                                                 
667   return fWrappedProcess->AtRestGetPhysicalInt    
668 }                                                 
669                                                   
670 G4VParticleChange* G4BiasingProcessInterface::    
671                                                   
672 {                                                 
673   return fWrappedProcess->AtRestDoIt(track, st    
674 }                                                 
675                                                   
676 G4bool G4BiasingProcessInterface::IsApplicable    
677 {                                                 
678   if ( fWrappedProcess != nullptr )               
679     return fWrappedProcess->IsApplicable(pd);     
680   else                                            
681     return true;                                  
682 }                                                 
683                                                   
684 void  G4BiasingProcessInterface::SetMasterProc    
685 {                                                 
686   // -- Master for this process:                  
687   G4VProcess::SetMasterProcess(masterP);          
688   // -- Master for wrapped process:               
689   if ( fWrappedProcess != nullptr )               
690   {                                               
691     const G4BiasingProcessInterface* thisWrapp    
692       = (const G4BiasingProcessInterface *)Get    
693     // -- paranoia check: (?)                     
694     G4VProcess* wrappedMaster = nullptr;          
695     wrappedMaster = thisWrapperMaster->GetWrap    
696     fWrappedProcess->SetMasterProcess( wrapped    
697   }                                               
698 }                                                 
699                                                   
700 void G4BiasingProcessInterface::                  
701 BuildPhysicsTable(const G4ParticleDefinition&     
702 {                                                 
703   // -- Sequential mode : called second (after    
704   // -- MT mode         : called second (after    
705   // --                   Corresponding proces    
706   // -- PreparePhysicsTable(...) has been call    
707   // -- so the first/last flags and G4BiasingP    
708   // -- been properly setup, fIamFirstGPIL is     
709   if ( fWrappedProcess != nullptr )               
710   {                                               
711     fWrappedProcess->BuildPhysicsTable(pd);       
712   }                                               
713                                                   
714   if ( fIamFirstGPIL )                            
715   {                                               
716     // -- Re-order vector of processes to matc    
717     // -- (made for fIamFirstGPIL, but importa    
718     ReorderBiasingVectorAsGPIL();                 
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    
740   SetUpFirstLastFlags();                          
741   if ( fWrappedProcess != nullptr )               
742   {                                               
743     fWrappedProcess->PreparePhysicsTable(pd);     
744   }                                               
745 }                                                 
746                                                   
747 G4bool  G4BiasingProcessInterface::               
748 StorePhysicsTable(const G4ParticleDefinition*     
749 {                                                 
750   if ( fWrappedProcess != nullptr )               
751     return fWrappedProcess->StorePhysicsTable(    
752   else                                            
753     return false;                                 
754 }                                                 
755                                                   
756 G4bool G4BiasingProcessInterface::                
757 RetrievePhysicsTable(const G4ParticleDefinitio    
758 {                                                 
759   if ( fWrappedProcess != nullptr )               
760     return fWrappedProcess->RetrievePhysicsTab    
761   else                                            
762     return false;                                 
763 }                                                 
764                                                   
765 void G4BiasingProcessInterface::SetProcessMana    
766 {                                                 
767   if ( fWrappedProcess != nullptr )               
768     fWrappedProcess->SetProcessManager(mgr);      
769   else                                            
770     G4VProcess::SetProcessManager(mgr);           
771                                                   
772   // -- initialize fSharedData pointer:           
773   if (G4BiasingProcessSharedData::fSharedDataM    
774    == G4BiasingProcessSharedData::fSharedDataM    
775   {                                               
776     fSharedData = new G4BiasingProcessSharedDa    
777     G4BiasingProcessSharedData::fSharedDataMap    
778   }                                               
779   else                                            
780   {                                               
781     fSharedData = G4BiasingProcessSharedData::    
782   }                                               
783   // -- augment list of co-operating processes    
784   fSharedData->fBiasingProcessInterfaces.push_    
785   fSharedData->fPublicBiasingProcessInterfaces    
786   if ( fIsPhysicsBasedBiasing )                   
787   {                                               
788     fSharedData->fPhysicsBiasingProcessInterfa    
789     fSharedData-> fPublicPhysicsBiasingProcess    
790   }                                               
791   else                                            
792   {                                               
793     fSharedData->fNonPhysicsBiasingProcessInte    
794     fSharedData->fPublicNonPhysicsBiasingProce    
795   }                                               
796   // -- remember process manager:                 
797   fProcessManager = mgr;                          
798 }                                                 
799                                                   
800 const G4ProcessManager* G4BiasingProcessInterf    
801 {                                                 
802   if ( fWrappedProcess != nullptr )               
803     return fWrappedProcess->GetProcessManager(    
804   else                                            
805     return G4VProcess::GetProcessManager();       
806 }                                                 
807                                                   
808 void G4BiasingProcessInterface::                  
809 BuildWorkerPhysicsTable(const G4ParticleDefini    
810 {                                                 
811   // -- Sequential mode : not called              
812   // -- MT mode         : called after Prepare    
813   // -- PrepareWorkerPhysicsTable(...) has bee    
814   // -- so the first/last flags and G4BiasingP    
815   // -- been properly setup, fIamFirstGPIL is     
816   if ( fWrappedProcess != nullptr )               
817   {                                               
818     fWrappedProcess->BuildWorkerPhysicsTable(p    
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     }                                             
837   }                                               
838 }                                                 
839                                                   
840 void G4BiasingProcessInterface::                  
841 PrepareWorkerPhysicsTable(const G4ParticleDefi    
842 {                                                 
843   // -- Sequential mode : not called              
844   // -- MT mode         : called first, before    
845   // -- Let process finding its first/last pos    
846   SetUpFirstLastFlags();                          
847                                                   
848   if ( fWrappedProcess != nullptr )               
849   {                                               
850     fWrappedProcess->PrepareWorkerPhysicsTable    
851   }                                               
852 }                                                 
853                                                   
854 void G4BiasingProcessInterface::ResetNumberOfI    
855 {                                                 
856   if ( fWrappedProcess != nullptr )               
857     fWrappedProcess->ResetNumberOfInteractionL    
858 }                                                 
859                                                   
860 G4bool G4BiasingProcessInterface::                
861 GetIsFirstPostStepGPILInterface( G4bool physOn    
862 {                                                 
863   G4int iPhys = ( physOnly ) ? 1 : 0;             
864   return fFirstLastFlags[IdxFirstLast( 1, 1, i    
865 }                                                 
866                                                   
867 G4bool  G4BiasingProcessInterface::               
868 GetIsLastPostStepGPILInterface( G4bool physOnl    
869 {                                                 
870   G4int iPhys = ( physOnly ) ? 1 : 0;             
871   return fFirstLastFlags[IdxFirstLast( 0, 1, i    
872 }                                                 
873                                                   
874 G4bool G4BiasingProcessInterface::                
875 GetIsFirstPostStepDoItInterface( G4bool physOn    
876 {                                                 
877   G4int iPhys = ( physOnly ) ? 1 : 0;             
878   return fFirstLastFlags[IdxFirstLast( 1, 0, i    
879 }                                                 
880                                                   
881 G4bool  G4BiasingProcessInterface::               
882 GetIsLastPostStepDoItInterface( G4bool physOnl    
883 {                                                 
884   G4int iPhys = ( physOnly ) ? 1 : 0;             
885   return fFirstLastFlags[IdxFirstLast( 0, 0, i    
886 }                                                 
887                                                   
888 G4bool G4BiasingProcessInterface::                
889 IsFirstPostStepGPILInterface(G4bool physOnly)     
890 {                                                 
891   G4bool isFirst = true;                          
892   const G4ProcessVector* pv = fProcessManager-    
893   G4int thisIdx(-1);                              
894   for ( auto i = 0; i < (G4int)pv->size(); ++i    
895   {                                               
896     if ( (*pv)(i) == this ) { thisIdx = i; bre    
897   }                                               
898   if ( thisIdx < 0 ) return false; // -- to ig    
899   for ( std::size_t i=0; i<(fSharedData->fBias    
900   {                                               
901     if ( (fSharedData->fBiasingProcessInterfac    
902     {                                             
903       G4int thatIdx(-1);                          
904       for (auto j = 0; j < (G4int)pv->size();     
905       {                                           
906         if ( (*pv)(j) == (fSharedData->fBiasin    
907         {                                         
908           thatIdx = j; break;                     
909         }                                         
910       }                                           
911       if ( thatIdx >= 0 ) // -- to ignore pure    
912       {                                           
913         if ( thisIdx >  thatIdx )                 
914         {                                         
915           isFirst = false;                        
916           break;                                  
917         }                                         
918       }                                           
919     }                                             
920   }                                               
921   return isFirst;                                 
922 }                                                 
923                                                   
924 G4bool G4BiasingProcessInterface::                
925 IsLastPostStepGPILInterface(G4bool physOnly) c    
926 {                                                 
927   G4bool isLast = true;                           
928   const G4ProcessVector* pv = fProcessManager-    
929   G4int thisIdx(-1);                              
930   for (auto i = 0; i < (G4int)pv->size(); ++i     
931   {                                               
932     if ( (*pv)(i) == this ) { thisIdx = i; bre    
933   }                                               
934   if ( thisIdx < 0 ) return false; // -- to ig    
935   for (std::size_t i=0; i<(fSharedData->fBiasi    
936   {                                               
937     if ( (fSharedData->fBiasingProcessInterfac    
938     {                                             
939       G4int thatIdx(-1);                          
940       for (auto j = 0; j < (G4int)pv->size();     
941       {                                           
942         if ( (*pv)(j) == (fSharedData->fBiasin    
943         {                                         
944           thatIdx = j; break;                     
945         }                                         
946       }                                           
947       if ( thatIdx >= 0 ) // -- to ignore pure    
948       {                                           
949         if ( thisIdx <  thatIdx )                 
950         {                                         
951           isLast = false;                         
952           break;                                  
953         }                                         
954       }                                           
955     }                                             
956   }                                               
957   return isLast;                                  
958 }                                                 
959                                                   
960 G4bool G4BiasingProcessInterface::                
961 IsFirstPostStepDoItInterface(G4bool physOnly)     
962 {                                                 
963   G4bool isFirst = true;                          
964   const G4ProcessVector* pv = fProcessManager-    
965   G4int thisIdx(-1);                              
966   for (auto i = 0; i < (G4int)pv->size(); ++i     
967   {                                               
968     if ( (*pv)(i) == this ) { thisIdx = i; bre    
969   }                                               
970   if ( thisIdx < 0 ) return false; // -- to ig    
971   for (std::size_t i=0; i<(fSharedData->fBiasi    
972   {                                               
973     if ( (fSharedData->fBiasingProcessInterfac    
974     {                                             
975       G4int thatIdx(-1);                          
976       for (auto j = 0; j < (G4int)pv->size();     
977       {                                           
978         if ( (*pv)(j) == (fSharedData->fBiasin    
979         {                                         
980           thatIdx = j; break;                     
981         }                                         
982       }                                           
983       if ( thatIdx >= 0 ) // -- to ignore pure    
984       {                                           
985         if ( thisIdx >  thatIdx )                 
986         {                                         
987           isFirst = false;                        
988           break;                                  
989         }                                         
990       }                                           
991     }                                             
992   }                                               
993   return isFirst;                                 
994 }                                                 
995                                                   
996 G4bool G4BiasingProcessInterface::                
997 IsLastPostStepDoItInterface(G4bool physOnly) c    
998 {                                                 
999   G4bool isLast = true;                           
1000   const G4ProcessVector* pv = fProcessManager    
1001   G4int thisIdx(-1);                             
1002   for (auto i = 0; i < (G4int)pv->size(); ++i    
1003   {                                              
1004     if ( (*pv)(i) == this ) { thisIdx = i; br    
1005   }                                              
1006   if ( thisIdx < 0 ) return false; // -- to i    
1007   for (std::size_t i=0; i<(fSharedData->fBias    
1008   {                                              
1009     if ( (fSharedData->fBiasingProcessInterfa    
1010     {                                            
1011       G4int thatIdx(-1);                         
1012       for (auto j = 0; j < (G4int)pv->size();    
1013       {                                          
1014         if ( (*pv)(j) == (fSharedData->fBiasi    
1015         {                                        
1016           thatIdx = j; break;                    
1017         }                                        
1018       }                                          
1019       if ( thatIdx >= 0 ) // -- to ignore pur    
1020       {                                          
1021         if ( thisIdx <  thatIdx )                
1022         {                                        
1023           isLast = false;                        
1024           break;                                 
1025         }                                        
1026       }                                          
1027     }                                            
1028   }                                              
1029   return isLast;                                 
1030 }                                                
1031                                                  
1032 void G4BiasingProcessInterface::SetUpFirstLas    
1033 {                                                
1034   for (G4int iPhys = 0; iPhys < 2; ++iPhys)      
1035   {                                              
1036     G4bool physOnly = ( iPhys == 1 );            
1037     fFirstLastFlags[IdxFirstLast( 1, 1, iPhys    
1038     fFirstLastFlags[IdxFirstLast( 0, 1, iPhys    
1039     fFirstLastFlags[IdxFirstLast( 1, 0, iPhys    
1040     fFirstLastFlags[IdxFirstLast( 0, 0, iPhys    
1041   }                                              
1042                                                  
1043   // -- for itself, for optimization:            
1044   fIamFirstGPIL =  GetIsFirstPostStepGPILInte    
1045 }                                                
1046                                                  
1047 void G4BiasingProcessInterface::ResetForUnbia    
1048 {                                                
1049   fOccurenceBiasingOperation  = nullptr;         
1050   fFinalStateBiasingOperation = nullptr;         
1051   fNonPhysicsBiasingOperation = nullptr;         
1052   fBiasingInteractionLaw      = nullptr;         
1053 }                                                
1054                                                  
1055 void G4BiasingProcessInterface::                 
1056 InvokeWrappedProcessPostStepGPIL( const G4Tra    
1057                                   G4double pr    
1058                                   G4ForceCond    
1059 {                                                
1060   G4double usedPreviousStepSize = previousSte    
1061   // -- if the physics process has been under    
1062   // -- we reset it, as we don't know if it w    
1063   // -- step. The pity is that PostStepGPIL a    
1064   // -- calculations are done both in the Pos    
1065   // -- are just interested in the calculatio    
1066   // -- as this forces to re-generated a rand    
1067   if ( fResetWrappedProcessInteractionLength     
1068   {                                              
1069     fResetWrappedProcessInteractionLength = f    
1070     fWrappedProcess->ResetNumberOfInteraction    
1071     // -- We set "previous step size" as 0.0,    
1072     usedPreviousStepSize = 0.0;                  
1073   }                                              
1074   // -- GPIL response:                           
1075   fWrappedProcessPostStepGPIL = fWrappedProce    
1076   fWrappedProcessForceCondition = *condition;    
1077   // -- and (inverse) cross-section:             
1078   fWrappedProcessInteractionLength = fWrapped    
1079 }                                                
1080                                                  
1081 void G4BiasingProcessInterface::ReorderBiasin    
1082 {                                                
1083   // -- re-order vector of processes to match    
1084   std::vector < G4BiasingProcessInterface* >     
1085   ( fSharedData -> fBiasingProcessInterfaces     
1086   ( fSharedData -> fPhysicsBiasingProcessInte    
1087   ( fSharedData -> fNonPhysicsBiasingProcessI    
1088   ( fSharedData -> fPublicBiasingProcessInter    
1089   ( fSharedData -> fPublicPhysicsBiasingProce    
1090   ( fSharedData -> fPublicNonPhysicsBiasingPr    
1091                                                  
1092   const G4ProcessVector* pv = fProcessManager    
1093   for (auto i = 0; i < (G4int)pv->size(); ++i    
1094   {                                              
1095     for (std::size_t j = 0; j < tmpProcess.si    
1096     {                                            
1097       if ( (*pv)(i) == tmpProcess[j] )           
1098       {                                          
1099         ( fSharedData->fBiasingProcessInterfa    
1100         ( fSharedData->fPublicBiasingProcessI    
1101         if ( tmpProcess[j] -> fIsPhysicsBased    
1102         {                                        
1103           ( fSharedData->fPhysicsBiasingProce    
1104           ( fSharedData->fPublicPhysicsBiasin    
1105         }                                        
1106         else                                     
1107         {                                        
1108           ( fSharedData -> fNonPhysicsBiasing    
1109           ( fSharedData -> fPublicNonPhysicsB    
1110         }                                        
1111         break;                                   
1112       }                                          
1113     }                                            
1114   }                                              
1115 }                                                
1116