| Geant4 Cross Reference |
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 // G4BOptnLeadingParticle
27 // -------------------------------------------
28
29 #include "G4BOptnLeadingParticle.hh"
30 #include "G4BiasingProcessInterface.hh"
31
32 #include <vector>
33 #include <map>
34
35
36 G4BOptnLeadingParticle::G4BOptnLeadingParticle
37 : G4VBiasingOperation( name )
38 {
39 }
40
41 G4BOptnLeadingParticle::~G4BOptnLeadingParticl
42 {
43 }
44
45 G4VParticleChange* G4BOptnLeadingParticle::
46 ApplyFinalStateBiasing( const G4BiasingProcess
47 const G4Track* track,
48 {
49 // -- collect wrapped process particle chang
50 auto wrappedProcessParticleChange = callingP
51
52 // -- does nothing in case the primary stays
53 if ( wrappedProcessParticleChange->GetNumber
54 if ( wrappedProcessParticleChange->GetTrackS
55 // -- ... else, collect the secondaries in a
56 std::vector < G4Track* > secondariesAndPrima
57 for ( auto i = 0 ; i < wrappedProcessParticl
58
59
60 // -- If case the primary survives, need to
61 // -- this is tricky, as this class does not
62 // -- caring about the touchables, etc.) So
63 // -- class with a dynamic cast. If we don't
64 // -- the untrimmed process final state.
65 // -- Note that this case does not happen of
66 // -- particles can be produced without kill
67 G4ParticleChange* castParticleChange ( nullp
68 G4Track* finalStatePrimary ( nullp
69 if ( ( wrappedProcessParticleChange->GetTrac
70 {
71 // fFakePrimaryTrack->CopyTrackInfo( *trac
72 // fFakeStep ->InitializeStep( fFak
73 // wrappedProcessParticleChange->UpdateSte
74 // fFakeStep->UpdateTrack();
75 castParticleChange = dynamic_cast< G4Parti
76 if ( castParticleChange == nullptr )
77 {
78 G4cout << " **** G4BOptnLeadingParticle:
79 return wrappedProcessParticleChange;
80 }
81 finalStatePrimary = new G4Track( *track );
82 finalStatePrimary->SetKineticEnergy ( cast
83 finalStatePrimary->SetWeight ( castParticl
84 finalStatePrimary->SetMomentumDirection( *
85 // -- [**] push the primary as the last tr
86 secondariesAndPrimary.push_back( finalStat
87 }
88
89 // -- Ensure the secondaries all have the pr
90 // ---- collect primary track weight, from u
91 G4double primaryWeight;
92 if ( finalStatePrimary ) primaryWeight = fin
93 else primaryWeight = track->GetWeight();
94 // ---- now set this same weight to all seco
95 for (auto i = 0; i < wrappedProcessParticleC
96 secondariesAndPrimary[ i ]->SetWeight( pri
97
98 // -- finds the leading particle, initialize
99 std::size_t leadingIDX = 0;
100 G4double leadingEnergy = -1;
101 std::map< G4Track*, G4bool > survivingMap;
102 for ( std::size_t idx = 0; idx < secondaries
103 {
104 survivingMap[ secondariesAndPrimary[idx] ]
105 if ( secondariesAndPrimary[idx]->GetKineti
106 {
107 leadingEnergy = secondariesAndPrimary[id
108 leadingIDX = idx;
109 }
110 }
111 survivingMap[ secondariesAndPrimary[leadingI
112
113 // -- now make track vectors of given types
114 std::map < G4int, std::vector< G4Track* > >
115 for ( std::size_t idx = 0; idx < secondaries
116 {
117 if ( idx == leadingIDX ) continue; // -- e
118 auto currentTrack = secondariesAndPrimary[
119 auto GROUP = std::abs( currentTrack
120 if ( currentTrack->GetDefinition()->GetBar
121 GROUP = -1000; // -- merge all baryons a
122
123 if ( typesAndTracks.find( GROUP ) == types
124 {
125 std::vector< G4Track* > v;
126 v.push_back( currentTrack );
127 typesAndTracks[ GROUP ] = std::move(v);
128 }
129 else
130 {
131 typesAndTracks[ GROUP ].push_back( curre
132 }
133 }
134 // -- and on these vectors, randomly select
135 // ---- randomly select one surviving track
136 // ---- for this surviving track, further ap
137 G4int nSecondaries = 0; // -- the number of
138 for ( auto& typeAndTrack : typesAndTracks )
139 {
140 std::size_t nTracks = (typeAndTrack.second
141 G4Track* keptTrack;
142 // -- select one track among ones in same
143 if ( nTracks > 1 )
144 {
145 auto keptTrackIDX = G4RandFlat::shootInt
146 keptTrack = (typeAndTrack.second)[keptTr
147 keptTrack->SetWeight( keptTrack->GetWeig
148 }
149 else
150 {
151 keptTrack = (typeAndTrack.second)[0];
152 }
153 // -- further apply a Russian Roulette on
154 G4bool keepTrack = false;
155 if ( fRussianRouletteKillingProbability >
156 {
157 if ( G4UniformRand() > fRussianRouletteK
158 {
159 keptTrack->SetWeight( keptTrack->GetWe
160 keepTrack = true;
161 }
162 }
163 else keepTrack = true;
164 if ( keepTrack )
165 {
166 survivingMap[ keptTrack ] = true;
167 if ( keptTrack != finalStatePrimary ) ++
168 }
169 }
170 // -- and if the leading is not the primary,
171 if ( secondariesAndPrimary[leadingIDX] != fi
172
173 // -- verify if the primary is still alive o
174 G4bool primarySurvived = false;
175 if ( finalStatePrimary )
176 primarySurvived = survivingMap[ finalState
177
178
179 // -- fill the trimmed particle change:
180 // ---- fill for the primary:
181 fParticleChange.Initialize(*track);
182 if ( primarySurvived )
183 {
184 fParticleChange.ProposeTrackStatus ( wrapp
185 fParticleChange.ProposeParentWeight ( fina
186 // -- take weight from copy of primary,
187 // random selection loop above
188 fParticleChange.ProposeEnergy ( finalState
189 fParticleChange.ProposeMomentumDirection (
190 }
191 else
192 {
193 fParticleChange.ProposeTrackStatus ( fStop
194 fParticleChange.ProposeParentWeight( 0.0 )
195 fParticleChange.ProposeEnergy ( 0.0 )
196 }
197 // -- fill for surviving secondaries:
198 fParticleChange.SetSecondaryWeightByProcess(
199 fParticleChange.SetNumberOfSecondaries(nSeco
200 // ---- note we loop up to on the number of
201 ////// G4cout << callingProcess->GetProcess
202 for ( auto idx = 0 ; idx < wrappedProcessPar
203 {
204 G4Track* secondary = secondariesAndPrimary
205 // ********************
206 //// if ( !survivingMap[ secondary ] ) G4c
207 ///// else G4cout << " ";
208 ///// G4cout << secondary->GetDefinition()
209 ///// if ( !survivingMap[ secondary ] ) G4
210 //// if ( secondary == secondariesAndPrim
211 // ******************
212 if ( survivingMap[ secondary ] )
213 fParticleChange.AddSecondary( secondary
214 else
215 delete secondary;
216 }
217 /// G4cout << G4endl;
218
219 // -- clean the wrapped process particle cha
220 wrappedProcessParticleChange->Clear();
221
222 if ( finalStatePrimary ) delete finalStatePr
223
224 // -- finally, returns the trimmed particle
225 return &fParticleChange;
226 }
227