G4BOptnLeadingParticle Class Reference

#include <G4BOptnLeadingParticle.hh>

Inheritance diagram for G4BOptnLeadingParticle:

Public Member Functions
	G4BOptnLeadingParticle (G4String name)
virtual	~G4BOptnLeadingParticle ()
virtual const G4VBiasingInteractionLaw *	ProvideOccurenceBiasingInteractionLaw (const G4BiasingProcessInterface *, G4ForceCondition &)
virtual G4VParticleChange *	ApplyFinalStateBiasing (const G4BiasingProcessInterface *, const G4Track *, const G4Step *, G4bool &)
virtual G4double	DistanceToApplyOperation (const G4Track *, G4double, G4ForceCondition *)
virtual G4VParticleChange *	GenerateBiasingFinalState (const G4Track *, const G4Step *)
void	SetFurtherKillingProbability (G4double p)
G4double	GetFurtherKillingProbability () const
Public Member Functions inherited from G4VBiasingOperation
	G4VBiasingOperation (G4String name)
virtual	~G4VBiasingOperation ()
virtual const G4VBiasingInteractionLaw *	ProvideOccurenceBiasingInteractionLaw (const G4BiasingProcessInterface *, G4ForceCondition &)=0
virtual G4double	ProposeAlongStepLimit (const G4BiasingProcessInterface *)
virtual G4GPILSelection	ProposeGPILSelection (const G4GPILSelection wrappedProcessSelection)
virtual void	AlongMoveBy (const G4BiasingProcessInterface *, const G4Step *, G4double)
virtual G4VParticleChange *	ApplyFinalStateBiasing (const G4BiasingProcessInterface *, const G4Track *, const G4Step *, G4bool &)=0
virtual G4double	DistanceToApplyOperation (const G4Track *, G4double, G4ForceCondition *)=0
virtual G4VParticleChange *	GenerateBiasingFinalState (const G4Track *, const G4Step *)=0
const G4String &	GetName () const
std::size_t	GetUniqueID () const

Detailed Description

Definition at line 54 of file G4BOptnLeadingParticle.hh.

Constructor & Destructor Documentation

G4BOptnLeadingParticle()

G4BOptnLeadingParticle::G4BOptnLeadingParticle

(

G4String

name

)

Definition at line 33 of file G4BOptnLeadingParticle.cc.

  : G4VBiasingOperation                ( name ),
    fRussianRouletteKillingProbability ( -1.0 )
{
}

~G4BOptnLeadingParticle()

G4BOptnLeadingParticle::~G4BOptnLeadingParticle ( )

virtual

Definition at line 39 of file G4BOptnLeadingParticle.cc.

{
}

Member Function Documentation

ApplyFinalStateBiasing()

G4VParticleChange * G4BOptnLeadingParticle::ApplyFinalStateBiasing ( const G4BiasingProcessInterface *  callingProcess, const G4Track *  track, const G4Step *  step, G4bool &    )

virtual

G4cout << G4endl;

Implements G4VBiasingOperation.

Definition at line 43 of file G4BOptnLeadingParticle.cc.

{
  // -- collect wrapped process particle change:
  auto wrappedProcessParticleChange = callingProcess->GetWrappedProcess()->PostStepDoIt(*track,*step);
 
  // -- does nothing in case the primary stays alone or in weird situation where all are killed...
  if ( wrappedProcessParticleChange->GetNumberOfSecondaries() == 0 )                        return wrappedProcessParticleChange;
  if ( wrappedProcessParticleChange->GetTrackStatus()         == fKillTrackAndSecondaries ) return wrappedProcessParticleChange;
  // -- ... else, collect the secondaries in a same vector (the primary is pushed in this vector, if surviving, later see [**]):
  std::vector < G4Track* > secondariesAndPrimary;
  for ( auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary.push_back(  wrappedProcessParticleChange->GetSecondary( i ) );
  
 
  // -- If case the primary survives, need to collect its new state. In the general case of the base class G4VParticleChange
  // -- this is tricky, as this class does not hold the primary changes (and we have to build a fake step and fake track
  // -- caring about the touchables, etc.) So we limit here to the G4ParticleChange case, checking the reality of this
  // -- class with a dynamic cast. If we don't have such an actual G4DynamicParticle object, we give up the biasing and return
  // -- the untrimmed process final state.
  // -- Note that this case does not happen often, as the technique is intended for inelastic processes. For case where several
  // -- particles can be produced without killing the primary, we have for example the electron-/positron-nuclear
  G4ParticleChange* castParticleChange ( nullptr );
  G4Track*          finalStatePrimary  ( nullptr );
  if ( ( wrappedProcessParticleChange->GetTrackStatus() != fStopAndKill ) )
    {
      // fFakePrimaryTrack->CopyTrackInfo( *track );
      // fFakeStep        ->InitializeStep( fFakePrimaryTrack );
      // wrappedProcessParticleChange->UpdateStepForPostStep( fFakeStep );
      // fFakeStep->UpdateTrack();
      castParticleChange = dynamic_cast< G4ParticleChange* >(  wrappedProcessParticleChange );
      if ( castParticleChange == nullptr )
    {
      G4cout << " **** G4BOptnLeadingParticle::ApplyFinalStateBiasing(...) : can not bias for " << callingProcess->GetProcessName() << ", this is just a warning." << G4endl;
      return wrappedProcessParticleChange;
    }
      finalStatePrimary = new G4Track( *track );
      finalStatePrimary->SetKineticEnergy    ( castParticleChange->GetEnergy()               );
      finalStatePrimary->SetWeight           ( castParticleChange->GetWeight()               );
      finalStatePrimary->SetMomentumDirection( *(castParticleChange->GetMomentumDirection()) );
      // -- [**] push the primary as the last track in the vector of tracks:
      secondariesAndPrimary.push_back( finalStatePrimary );
    }
  
  // -- Ensure the secondaries all have the primary weight:
  // ---- collect primary track weight, from updated by process if alive, or from original copy if died:
  G4double                 primaryWeight;
  if ( finalStatePrimary ) primaryWeight = finalStatePrimary->GetWeight();
  else                     primaryWeight = track            ->GetWeight();
  // ---- now set this same weight to all secondaries:
  for ( auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary[ i ]->SetWeight( primaryWeight );
  
 
  // -- finds the leading particle, initialize a map of surviving tracks, tag as surviving the leading track:
  size_t   leadingIDX    =  0;
  G4double leadingEnergy = -1;
  std::map< G4Track*, G4bool > survivingMap;
  for ( size_t idx = 0; idx < secondariesAndPrimary.size(); idx++ )
    {
      survivingMap[ secondariesAndPrimary[idx] ] = false;
      if ( secondariesAndPrimary[idx]->GetKineticEnergy() > leadingEnergy )
    {
      leadingEnergy = secondariesAndPrimary[idx]->GetKineticEnergy();
      leadingIDX    = idx;
    }
    }
  survivingMap[ secondariesAndPrimary[leadingIDX] ] = true; // -- tag as surviving the leading particle
  
  
  // -- now make track vectors of given types ( choose type = abs(PDG) ), excluding the leading particle:
  std::map < G4int, std::vector< G4Track* > > typesAndTracks;
  for ( size_t idx = 0; idx < secondariesAndPrimary.size(); idx++ )
    {
      if ( idx == leadingIDX ) continue; // -- excludes the leading particle
      auto currentTrack = secondariesAndPrimary[idx];
      auto GROUP        = std::abs( currentTrack->GetDefinition()->GetPDGEncoding() ); // -- merge particles and anti-particles in the same category  -- §§ this might be proposed as an option in future
      if ( currentTrack->GetDefinition()->GetBaryonNumber() >= 2 ) GROUP = -1000; // -- merge all baryons above proton/neutron in one same group -- §§ might be proposed as an option too
      
      if ( typesAndTracks.find( GROUP ) == typesAndTracks.end() )
    {
      std::vector< G4Track* > v;
      v.push_back( currentTrack );
      typesAndTracks[ GROUP ] = v;
    }
      else
    {
      typesAndTracks[ GROUP ].push_back( currentTrack );
    }
    }
  // -- and on these vectors, randomly select the surviving particles:
  // ---- randomly select one surviving track per species
  // ---- for this surviving track, further apply a Russian roulette
  G4int nSecondaries = 0; // -- the number of secondaries to be returned
  for ( auto typeAndTrack : typesAndTracks )
    {
      size_t nTracks = (typeAndTrack.second).size();
      G4Track* keptTrack;
      // -- select one track among ones in same species:
      if ( nTracks > 1 )
    {
      auto keptTrackIDX = G4RandFlat::shootInt( nTracks );
      keptTrack = (typeAndTrack.second)[keptTrackIDX];
      keptTrack->SetWeight( keptTrack->GetWeight() * nTracks );
    }
      else
    {
      keptTrack = (typeAndTrack.second)[0];
    }
      // -- further apply a Russian Roulette on it:
      G4bool keepTrack = false;
      if ( fRussianRouletteKillingProbability > 0.0 )
    {
      if ( G4UniformRand() > fRussianRouletteKillingProbability )
        {
          keptTrack->SetWeight( keptTrack->GetWeight() / (1. - fRussianRouletteKillingProbability) );
          keepTrack = true;
        }
    }
      else keepTrack = true;
      if ( keepTrack )
    {
      survivingMap[ keptTrack ] = true;
      if ( keptTrack != finalStatePrimary ) nSecondaries++;
    }
    }
  // -- and if the leading is not the primary, we have to count it in nSecondaries:
  if ( secondariesAndPrimary[leadingIDX] != finalStatePrimary ) nSecondaries++;
  
  // -- verify if the primary is still alive or not after above selection:
  G4bool primarySurvived = false;
  if ( finalStatePrimary ) primarySurvived = survivingMap[ finalStatePrimary ];
    
  
  // -- fill the trimmed particle change:
  // ---- fill for the primary:
  fParticleChange.Initialize(*track);
  if ( primarySurvived )
    {
      fParticleChange.ProposeTrackStatus       ( wrappedProcessParticleChange->GetTrackStatus()   );
      fParticleChange.ProposeParentWeight      ( finalStatePrimary->GetWeight()                   ); // -- take weight from copy of primary, this one being updated in the random selection loop above
      fParticleChange.ProposeEnergy            ( finalStatePrimary->GetKineticEnergy()            );
      fParticleChange.ProposeMomentumDirection ( finalStatePrimary->GetMomentumDirection()        );
    }
  else
    {
      fParticleChange.ProposeTrackStatus ( fStopAndKill );
      fParticleChange.ProposeParentWeight( 0.0 );
      fParticleChange.ProposeEnergy      ( 0.0 );
    }
  // -- fill for surviving secondaries:
  fParticleChange.SetSecondaryWeightByProcess(true);
  fParticleChange.SetNumberOfSecondaries(nSecondaries);
  // ---- note we loop up to on the number of secondaries, which excludes the primary, last in secondariesAndPrimary vector:
  //////  G4cout << callingProcess->GetProcessName() << " :";
  for ( auto idx = 0 ; idx < wrappedProcessParticleChange->GetNumberOfSecondaries() ; idx++ )
    {
      G4Track* secondary = secondariesAndPrimary[idx];
      // ********************
      //// if ( !survivingMap[ secondary ] ) G4cout << " [";
      ///// else G4cout << " ";
      ///// G4cout << secondary->GetDefinition()->GetParticleName() << " " << secondary->GetKineticEnergy();
      ///// if ( !survivingMap[ secondary ] ) G4cout << "]";
      //// if ( secondary ==  secondariesAndPrimary[leadingIDX] ) G4cout << " ***";
      // ******************
      if ( survivingMap[ secondary ] )  fParticleChange.AddSecondary( secondary );
      else                              delete secondary;
    }
  ///  G4cout << G4endl;
  
  // -- clean the wrapped process particle change:
  wrappedProcessParticleChange->Clear();
 
  if ( finalStatePrimary ) delete finalStatePrimary;  
 
  // -- finally, returns the trimmed particle change:
  return &fParticleChange;
  
}

DistanceToApplyOperation()

virtual G4double G4BOptnLeadingParticle::DistanceToApplyOperation ( const G4Track *  , G4double  , G4ForceCondition *    )

inlinevirtual

Implements G4VBiasingOperation.

Definition at line 72 of file G4BOptnLeadingParticle.hh.

                                                                                               {return 0;}

GenerateBiasingFinalState()

virtual G4VParticleChange * G4BOptnLeadingParticle::GenerateBiasingFinalState ( const G4Track *  , const G4Step *    )

inlinevirtual

Implements G4VBiasingOperation.

Definition at line 75 of file G4BOptnLeadingParticle.hh.

                                                                                               {return nullptr;}

GetFurtherKillingProbability()

G4double G4BOptnLeadingParticle::GetFurtherKillingProbability ( ) const

inline

Definition at line 86 of file G4BOptnLeadingParticle.hh.

{ return fRussianRouletteKillingProbability; }

ProvideOccurenceBiasingInteractionLaw()

virtual const G4VBiasingInteractionLaw * G4BOptnLeadingParticle::ProvideOccurenceBiasingInteractionLaw ( const G4BiasingProcessInterface *  , G4ForceCondition &    )

inlinevirtual

Implements G4VBiasingOperation.

Definition at line 65 of file G4BOptnLeadingParticle.hh.

{return nullptr;}

SetFurtherKillingProbability()

void G4BOptnLeadingParticle::SetFurtherKillingProbability ( G4double  p )

inline

Definition at line 85 of file G4BOptnLeadingParticle.hh.

{ fRussianRouletteKillingProbability = p;    }   // -- if p <= 0.0 the killing is ignored.

---

The documentation for this class was generated from the following files:

• G4BOptnLeadingParticle.hh
• G4BOptnLeadingParticle.cc