| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4BOptnLeadingParticle <<
27 // ------------------------------------------- <<
28 <<
29 #include "G4BOptnLeadingParticle.hh" 26 #include "G4BOptnLeadingParticle.hh"
30 #include "G4BiasingProcessInterface.hh" 27 #include "G4BiasingProcessInterface.hh"
31 28
32 #include <vector> 29 #include <vector>
33 #include <map> 30 #include <map>
34 31
35 32
36 G4BOptnLeadingParticle::G4BOptnLeadingParticle << 33 G4BOptnLeadingParticle::G4BOptnLeadingParticle(G4String name)
37 : G4VBiasingOperation( name ) << 34 : G4VBiasingOperation ( name ),
>> 35 fRussianRouletteKillingProbability ( -1.0 )
38 { 36 {
39 } 37 }
40 38
41 G4BOptnLeadingParticle::~G4BOptnLeadingParticl 39 G4BOptnLeadingParticle::~G4BOptnLeadingParticle()
42 { 40 {
43 } 41 }
44 42
45 G4VParticleChange* G4BOptnLeadingParticle:: << 43 G4VParticleChange* G4BOptnLeadingParticle::ApplyFinalStateBiasing( const G4BiasingProcessInterface* callingProcess,
46 ApplyFinalStateBiasing( const G4BiasingProcess << 44 const G4Track* track,
47 const G4Track* track, << 45 const G4Step* step,
>> 46 G4bool& )
48 { 47 {
49 // -- collect wrapped process particle chang 48 // -- collect wrapped process particle change:
50 auto wrappedProcessParticleChange = callingP 49 auto wrappedProcessParticleChange = callingProcess->GetWrappedProcess()->PostStepDoIt(*track,*step);
51 50
52 // -- does nothing in case the primary stays 51 // -- does nothing in case the primary stays alone or in weird situation where all are killed...
53 if ( wrappedProcessParticleChange->GetNumber 52 if ( wrappedProcessParticleChange->GetNumberOfSecondaries() == 0 ) return wrappedProcessParticleChange;
54 if ( wrappedProcessParticleChange->GetTrackS 53 if ( wrappedProcessParticleChange->GetTrackStatus() == fKillTrackAndSecondaries ) return wrappedProcessParticleChange;
55 // -- ... else, collect the secondaries in a 54 // -- ... else, collect the secondaries in a same vector (the primary is pushed in this vector, if surviving, later see [**]):
56 std::vector < G4Track* > secondariesAndPrima 55 std::vector < G4Track* > secondariesAndPrimary;
57 for ( auto i = 0 ; i < wrappedProcessParticl 56 for ( auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary.push_back( wrappedProcessParticleChange->GetSecondary( i ) );
58 57
59 58
60 // -- If case the primary survives, need to 59 // -- If case the primary survives, need to collect its new state. In the general case of the base class G4VParticleChange
61 // -- this is tricky, as this class does not 60 // -- this is tricky, as this class does not hold the primary changes (and we have to build a fake step and fake track
62 // -- caring about the touchables, etc.) So 61 // -- caring about the touchables, etc.) So we limit here to the G4ParticleChange case, checking the reality of this
63 // -- class with a dynamic cast. If we don't 62 // -- class with a dynamic cast. If we don't have such an actual G4DynamicParticle object, we give up the biasing and return
64 // -- the untrimmed process final state. 63 // -- the untrimmed process final state.
65 // -- Note that this case does not happen of 64 // -- Note that this case does not happen often, as the technique is intended for inelastic processes. For case where several
66 // -- particles can be produced without kill 65 // -- particles can be produced without killing the primary, we have for example the electron-/positron-nuclear
67 G4ParticleChange* castParticleChange ( nullp 66 G4ParticleChange* castParticleChange ( nullptr );
68 G4Track* finalStatePrimary ( nullp 67 G4Track* finalStatePrimary ( nullptr );
69 if ( ( wrappedProcessParticleChange->GetTrac 68 if ( ( wrappedProcessParticleChange->GetTrackStatus() != fStopAndKill ) )
70 { << 69 {
71 // fFakePrimaryTrack->CopyTrackInfo( *trac << 70 // fFakePrimaryTrack->CopyTrackInfo( *track );
72 // fFakeStep ->InitializeStep( fFak << 71 // fFakeStep ->InitializeStep( fFakePrimaryTrack );
73 // wrappedProcessParticleChange->UpdateSte << 72 // wrappedProcessParticleChange->UpdateStepForPostStep( fFakeStep );
74 // fFakeStep->UpdateTrack(); << 73 // fFakeStep->UpdateTrack();
75 castParticleChange = dynamic_cast< G4Parti << 74 castParticleChange = dynamic_cast< G4ParticleChange* >( wrappedProcessParticleChange );
76 if ( castParticleChange == nullptr ) << 75 if ( castParticleChange == nullptr )
77 { << 76 {
78 G4cout << " **** G4BOptnLeadingParticle: << 77 G4cout << " **** G4BOptnLeadingParticle::ApplyFinalStateBiasing(...) : can not bias for " << callingProcess->GetProcessName() << ", this is just a warning." << G4endl;
79 return wrappedProcessParticleChange; << 78 return wrappedProcessParticleChange;
80 } << 79 }
81 finalStatePrimary = new G4Track( *track ); << 80 finalStatePrimary = new G4Track( *track );
82 finalStatePrimary->SetKineticEnergy ( cast << 81 finalStatePrimary->SetKineticEnergy ( castParticleChange->GetEnergy() );
83 finalStatePrimary->SetWeight ( castParticl << 82 finalStatePrimary->SetWeight ( castParticleChange->GetWeight() );
84 finalStatePrimary->SetMomentumDirection( * << 83 finalStatePrimary->SetMomentumDirection( *(castParticleChange->GetMomentumDirection()) );
85 // -- [**] push the primary as the last tr << 84 // -- [**] push the primary as the last track in the vector of tracks:
86 secondariesAndPrimary.push_back( finalStat << 85 secondariesAndPrimary.push_back( finalStatePrimary );
87 } << 86 }
88 87
89 // -- Ensure the secondaries all have the pr 88 // -- Ensure the secondaries all have the primary weight:
90 // ---- collect primary track weight, from u 89 // ---- collect primary track weight, from updated by process if alive, or from original copy if died:
91 G4double primaryWeight; << 90 G4double primaryWeight;
92 if ( finalStatePrimary ) primaryWeight = fin 91 if ( finalStatePrimary ) primaryWeight = finalStatePrimary->GetWeight();
93 else primaryWeight = track->GetWeight(); << 92 else primaryWeight = track ->GetWeight();
94 // ---- now set this same weight to all seco 93 // ---- now set this same weight to all secondaries:
95 for (auto i = 0; i < wrappedProcessParticleC << 94 for ( auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary[ i ]->SetWeight( primaryWeight );
96 secondariesAndPrimary[ i ]->SetWeight( pri << 95
97 96
98 // -- finds the leading particle, initialize 97 // -- finds the leading particle, initialize a map of surviving tracks, tag as surviving the leading track:
99 std::size_t leadingIDX = 0; << 98 size_t leadingIDX = 0;
100 G4double leadingEnergy = -1; 99 G4double leadingEnergy = -1;
101 std::map< G4Track*, G4bool > survivingMap; 100 std::map< G4Track*, G4bool > survivingMap;
102 for ( std::size_t idx = 0; idx < secondaries << 101 for ( size_t idx = 0; idx < secondariesAndPrimary.size(); idx++ )
103 { <<
104 survivingMap[ secondariesAndPrimary[idx] ] <<
105 if ( secondariesAndPrimary[idx]->GetKineti <<
106 { 102 {
107 leadingEnergy = secondariesAndPrimary[id << 103 survivingMap[ secondariesAndPrimary[idx] ] = false;
108 leadingIDX = idx; << 104 if ( secondariesAndPrimary[idx]->GetKineticEnergy() > leadingEnergy )
>> 105 {
>> 106 leadingEnergy = secondariesAndPrimary[idx]->GetKineticEnergy();
>> 107 leadingIDX = idx;
>> 108 }
109 } 109 }
110 } <<
111 survivingMap[ secondariesAndPrimary[leadingI 110 survivingMap[ secondariesAndPrimary[leadingIDX] ] = true; // -- tag as surviving the leading particle
112 << 111
>> 112
113 // -- now make track vectors of given types 113 // -- now make track vectors of given types ( choose type = abs(PDG) ), excluding the leading particle:
114 std::map < G4int, std::vector< G4Track* > > 114 std::map < G4int, std::vector< G4Track* > > typesAndTracks;
115 for ( std::size_t idx = 0; idx < secondaries << 115 for ( size_t idx = 0; idx < secondariesAndPrimary.size(); idx++ )
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 { 116 {
125 std::vector< G4Track* > v; << 117 if ( idx == leadingIDX ) continue; // -- excludes the leading particle
126 v.push_back( currentTrack ); << 118 auto currentTrack = secondariesAndPrimary[idx];
127 typesAndTracks[ GROUP ] = std::move(v); << 119 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
>> 120 if ( currentTrack->GetDefinition()->GetBaryonNumber() >= 2 ) GROUP = -1000; // -- merge all baryons above proton/neutron in one same group -- §§ might be proposed as an option too
>> 121
>> 122 if ( typesAndTracks.find( GROUP ) == typesAndTracks.end() )
>> 123 {
>> 124 std::vector< G4Track* > v;
>> 125 v.push_back( currentTrack );
>> 126 typesAndTracks[ GROUP ] = v;
>> 127 }
>> 128 else
>> 129 {
>> 130 typesAndTracks[ GROUP ].push_back( currentTrack );
>> 131 }
128 } 132 }
129 else <<
130 { <<
131 typesAndTracks[ GROUP ].push_back( curre <<
132 } <<
133 } <<
134 // -- and on these vectors, randomly select 133 // -- and on these vectors, randomly select the surviving particles:
135 // ---- randomly select one surviving track 134 // ---- randomly select one surviving track per species
136 // ---- for this surviving track, further ap 135 // ---- for this surviving track, further apply a Russian roulette
137 G4int nSecondaries = 0; // -- the number of 136 G4int nSecondaries = 0; // -- the number of secondaries to be returned
138 for ( auto& typeAndTrack : typesAndTracks ) << 137 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 { 138 {
166 survivingMap[ keptTrack ] = true; << 139 size_t nTracks = (typeAndTrack.second).size();
167 if ( keptTrack != finalStatePrimary ) ++ << 140 G4Track* keptTrack;
>> 141 // -- select one track among ones in same species:
>> 142 if ( nTracks > 1 )
>> 143 {
>> 144 auto keptTrackIDX = G4RandFlat::shootInt( nTracks );
>> 145 keptTrack = (typeAndTrack.second)[keptTrackIDX];
>> 146 keptTrack->SetWeight( keptTrack->GetWeight() * nTracks );
>> 147 }
>> 148 else
>> 149 {
>> 150 keptTrack = (typeAndTrack.second)[0];
>> 151 }
>> 152 // -- further apply a Russian Roulette on it:
>> 153 G4bool keepTrack = false;
>> 154 if ( fRussianRouletteKillingProbability > 0.0 )
>> 155 {
>> 156 if ( G4UniformRand() > fRussianRouletteKillingProbability )
>> 157 {
>> 158 keptTrack->SetWeight( keptTrack->GetWeight() / (1. - fRussianRouletteKillingProbability) );
>> 159 keepTrack = true;
>> 160 }
>> 161 }
>> 162 else keepTrack = true;
>> 163 if ( keepTrack )
>> 164 {
>> 165 survivingMap[ keptTrack ] = true;
>> 166 if ( keptTrack != finalStatePrimary ) nSecondaries++;
>> 167 }
168 } 168 }
169 } <<
170 // -- and if the leading is not the primary, 169 // -- and if the leading is not the primary, we have to count it in nSecondaries:
171 if ( secondariesAndPrimary[leadingIDX] != fi << 170 if ( secondariesAndPrimary[leadingIDX] != finalStatePrimary ) nSecondaries++;
172 171
173 // -- verify if the primary is still alive o 172 // -- verify if the primary is still alive or not after above selection:
174 G4bool primarySurvived = false; 173 G4bool primarySurvived = false;
175 if ( finalStatePrimary ) << 174 if ( finalStatePrimary ) primarySurvived = survivingMap[ finalStatePrimary ];
176 primarySurvived = survivingMap[ finalState <<
177 175
178 176
179 // -- fill the trimmed particle change: 177 // -- fill the trimmed particle change:
180 // ---- fill for the primary: 178 // ---- fill for the primary:
181 fParticleChange.Initialize(*track); 179 fParticleChange.Initialize(*track);
182 if ( primarySurvived ) 180 if ( primarySurvived )
183 { << 181 {
184 fParticleChange.ProposeTrackStatus ( wrapp << 182 fParticleChange.ProposeTrackStatus ( wrappedProcessParticleChange->GetTrackStatus() );
185 fParticleChange.ProposeParentWeight ( fina << 183 fParticleChange.ProposeParentWeight ( finalStatePrimary->GetWeight() ); // -- take weight from copy of primary, this one being updated in the random selection loop above
186 // -- take weight from copy of primary, << 184 fParticleChange.ProposeEnergy ( finalStatePrimary->GetKineticEnergy() );
187 // random selection loop above << 185 fParticleChange.ProposeMomentumDirection ( finalStatePrimary->GetMomentumDirection() );
188 fParticleChange.ProposeEnergy ( finalState << 186 }
189 fParticleChange.ProposeMomentumDirection ( <<
190 } <<
191 else 187 else
192 { << 188 {
193 fParticleChange.ProposeTrackStatus ( fStop << 189 fParticleChange.ProposeTrackStatus ( fStopAndKill );
194 fParticleChange.ProposeParentWeight( 0.0 ) << 190 fParticleChange.ProposeParentWeight( 0.0 );
195 fParticleChange.ProposeEnergy ( 0.0 ) << 191 fParticleChange.ProposeEnergy ( 0.0 );
196 } 192 }
197 // -- fill for surviving secondaries: 193 // -- fill for surviving secondaries:
198 fParticleChange.SetSecondaryWeightByProcess( 194 fParticleChange.SetSecondaryWeightByProcess(true);
199 fParticleChange.SetNumberOfSecondaries(nSeco 195 fParticleChange.SetNumberOfSecondaries(nSecondaries);
200 // ---- note we loop up to on the number of 196 // ---- note we loop up to on the number of secondaries, which excludes the primary, last in secondariesAndPrimary vector:
201 ////// G4cout << callingProcess->GetProcess 197 ////// G4cout << callingProcess->GetProcessName() << " :";
202 for ( auto idx = 0 ; idx < wrappedProcessPar << 198 for ( auto idx = 0 ; idx < wrappedProcessParticleChange->GetNumberOfSecondaries() ; idx++ )
203 { << 199 {
204 G4Track* secondary = secondariesAndPrimary << 200 G4Track* secondary = secondariesAndPrimary[idx];
205 // ******************** << 201 // ********************
206 //// if ( !survivingMap[ secondary ] ) G4c << 202 //// if ( !survivingMap[ secondary ] ) G4cout << " [";
207 ///// else G4cout << " "; << 203 ///// else G4cout << " ";
208 ///// G4cout << secondary->GetDefinition() << 204 ///// G4cout << secondary->GetDefinition()->GetParticleName() << " " << secondary->GetKineticEnergy();
209 ///// if ( !survivingMap[ secondary ] ) G4 << 205 ///// if ( !survivingMap[ secondary ] ) G4cout << "]";
210 //// if ( secondary == secondariesAndPrim << 206 //// if ( secondary == secondariesAndPrimary[leadingIDX] ) G4cout << " ***";
211 // ****************** << 207 // ******************
212 if ( survivingMap[ secondary ] ) << 208 if ( survivingMap[ secondary ] ) fParticleChange.AddSecondary( secondary );
213 fParticleChange.AddSecondary( secondary << 209 else delete secondary;
214 else << 210 }
215 delete secondary; <<
216 } <<
217 /// G4cout << G4endl; 211 /// G4cout << G4endl;
218 212
219 // -- clean the wrapped process particle cha 213 // -- clean the wrapped process particle change:
220 wrappedProcessParticleChange->Clear(); 214 wrappedProcessParticleChange->Clear();
221 215
222 if ( finalStatePrimary ) delete finalStatePr 216 if ( finalStatePrimary ) delete finalStatePrimary;
223 217
224 // -- finally, returns the trimmed particle 218 // -- finally, returns the trimmed particle change:
225 return &fParticleChange; 219 return &fParticleChange;
>> 220
226 } 221 }
227 222