35 fRussianRouletteKillingProbability ( -1.0 )
52 if ( wrappedProcessParticleChange->GetNumberOfSecondaries() == 0 )
return wrappedProcessParticleChange;
53 if ( wrappedProcessParticleChange->GetTrackStatus() ==
fKillTrackAndSecondaries )
return wrappedProcessParticleChange;
55 std::vector < G4Track* > secondariesAndPrimary;
56 for (
auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary.push_back( wrappedProcessParticleChange->GetSecondary( i ) );
67 G4Track* finalStatePrimary (
nullptr );
68 if ( ( wrappedProcessParticleChange->GetTrackStatus() !=
fStopAndKill ) )
74 castParticleChange =
dynamic_cast< G4ParticleChange*
>( wrappedProcessParticleChange );
75 if ( castParticleChange ==
nullptr )
77 G4cout <<
" **** G4BOptnLeadingParticle::ApplyFinalStateBiasing(...) : can not bias for " << callingProcess->
GetProcessName() <<
", this is just a warning." <<
G4endl;
78 return wrappedProcessParticleChange;
80 finalStatePrimary =
new G4Track( *track );
85 secondariesAndPrimary.push_back( finalStatePrimary );
91 if ( finalStatePrimary ) primaryWeight = finalStatePrimary->
GetWeight();
94 for (
auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary[ i ]->SetWeight( primaryWeight );
98 size_t leadingIDX = 0;
100 std::map< G4Track*, G4bool > survivingMap;
101 for (
size_t idx = 0;
idx < secondariesAndPrimary.size();
idx++ )
103 survivingMap[ secondariesAndPrimary[
idx] ] =
false;
104 if ( secondariesAndPrimary[
idx]->GetKineticEnergy() > leadingEnergy )
106 leadingEnergy = secondariesAndPrimary[
idx]->GetKineticEnergy();
110 survivingMap[ secondariesAndPrimary[leadingIDX] ] =
true;
114 std::map < G4int, std::vector< G4Track* > > typesAndTracks;
115 for (
size_t idx = 0;
idx < secondariesAndPrimary.size();
idx++ )
117 if (
idx == leadingIDX )
continue;
118 auto currentTrack = secondariesAndPrimary[
idx];
119 auto GROUP =
std::abs( currentTrack->GetDefinition()->GetPDGEncoding() );
120 if ( currentTrack->GetDefinition()->GetBaryonNumber() >= 2 ) GROUP = -1000;
122 if ( typesAndTracks.find( GROUP ) == typesAndTracks.end() )
124 std::vector< G4Track* >
v;
125 v.push_back( currentTrack );
126 typesAndTracks[ GROUP ] =
v;
130 typesAndTracks[ GROUP ].push_back( currentTrack );
136 G4int nSecondaries = 0;
137 for (
auto typeAndTrack : typesAndTracks )
139 size_t nTracks = (typeAndTrack.second).size();
144 auto keptTrackIDX = G4RandFlat::shootInt( nTracks );
145 keptTrack = (typeAndTrack.second)[keptTrackIDX];
150 keptTrack = (typeAndTrack.second)[0];
162 else keepTrack =
true;
165 survivingMap[ keptTrack ] =
true;
166 if ( keptTrack != finalStatePrimary ) nSecondaries++;
170 if ( secondariesAndPrimary[leadingIDX] != finalStatePrimary ) nSecondaries++;
173 G4bool primarySurvived =
false;
174 if ( finalStatePrimary ) primarySurvived = survivingMap[ finalStatePrimary ];
180 if ( primarySurvived )
198 for (
auto idx = 0 ;
idx < wrappedProcessParticleChange->GetNumberOfSecondaries() ;
idx++ )
200 G4Track* secondary = secondariesAndPrimary[
idx];
209 else delete secondary;
214 wrappedProcessParticleChange->Clear();
216 if ( finalStatePrimary )
delete finalStatePrimary;