Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 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(const G4String& name) 37 : G4VBiasingOperation( name ) 38 { 39 } 40 41 G4BOptnLeadingParticle::~G4BOptnLeadingParticle() 42 { 43 } 44 45 G4VParticleChange* G4BOptnLeadingParticle:: 46 ApplyFinalStateBiasing( const G4BiasingProcessInterface* callingProcess, 47 const G4Track* track, const G4Step* step, G4bool& ) 48 { 49 // -- collect wrapped process particle change: 50 auto wrappedProcessParticleChange = callingProcess->GetWrappedProcess()->PostStepDoIt(*track,*step); 51 52 // -- does nothing in case the primary stays alone or in weird situation where all are killed... 53 if ( wrappedProcessParticleChange->GetNumberOfSecondaries() == 0 ) return wrappedProcessParticleChange; 54 if ( wrappedProcessParticleChange->GetTrackStatus() == fKillTrackAndSecondaries ) return wrappedProcessParticleChange; 55 // -- ... else, collect the secondaries in a same vector (the primary is pushed in this vector, if surviving, later see [**]): 56 std::vector < G4Track* > secondariesAndPrimary; 57 for ( auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary.push_back( wrappedProcessParticleChange->GetSecondary( i ) ); 58 59 60 // -- 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 hold the primary changes (and we have to build a fake step and fake track 62 // -- 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 have such an actual G4DynamicParticle object, we give up the biasing and return 64 // -- the untrimmed process final state. 65 // -- 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 killing the primary, we have for example the electron-/positron-nuclear 67 G4ParticleChange* castParticleChange ( nullptr ); 68 G4Track* finalStatePrimary ( nullptr ); 69 if ( ( wrappedProcessParticleChange->GetTrackStatus() != fStopAndKill ) ) 70 { 71 // fFakePrimaryTrack->CopyTrackInfo( *track ); 72 // fFakeStep ->InitializeStep( fFakePrimaryTrack ); 73 // wrappedProcessParticleChange->UpdateStepForPostStep( fFakeStep ); 74 // fFakeStep->UpdateTrack(); 75 castParticleChange = dynamic_cast< G4ParticleChange* >( wrappedProcessParticleChange ); 76 if ( castParticleChange == nullptr ) 77 { 78 G4cout << " **** G4BOptnLeadingParticle::ApplyFinalStateBiasing(...) : can not bias for " << callingProcess->GetProcessName() << ", this is just a warning." << G4endl; 79 return wrappedProcessParticleChange; 80 } 81 finalStatePrimary = new G4Track( *track ); 82 finalStatePrimary->SetKineticEnergy ( castParticleChange->GetEnergy() ); 83 finalStatePrimary->SetWeight ( castParticleChange->GetWeight() ); 84 finalStatePrimary->SetMomentumDirection( *(castParticleChange->GetMomentumDirection()) ); 85 // -- [**] push the primary as the last track in the vector of tracks: 86 secondariesAndPrimary.push_back( finalStatePrimary ); 87 } 88 89 // -- Ensure the secondaries all have the primary weight: 90 // ---- collect primary track weight, from updated by process if alive, or from original copy if died: 91 G4double primaryWeight; 92 if ( finalStatePrimary ) primaryWeight = finalStatePrimary->GetWeight(); 93 else primaryWeight = track->GetWeight(); 94 // ---- now set this same weight to all secondaries: 95 for (auto i = 0; i < wrappedProcessParticleChange->GetNumberOfSecondaries(); ++i) 96 secondariesAndPrimary[ i ]->SetWeight( primaryWeight ); 97 98 // -- finds the leading particle, initialize a map of surviving tracks, tag as surviving the leading track: 99 std::size_t leadingIDX = 0; 100 G4double leadingEnergy = -1; 101 std::map< G4Track*, G4bool > survivingMap; 102 for ( std::size_t idx = 0; idx < secondariesAndPrimary.size(); ++idx ) 103 { 104 survivingMap[ secondariesAndPrimary[idx] ] = false; 105 if ( secondariesAndPrimary[idx]->GetKineticEnergy() > leadingEnergy ) 106 { 107 leadingEnergy = secondariesAndPrimary[idx]->GetKineticEnergy(); 108 leadingIDX = idx; 109 } 110 } 111 survivingMap[ secondariesAndPrimary[leadingIDX] ] = true; // -- tag as surviving the leading particle 112 113 // -- now make track vectors of given types ( choose type = abs(PDG) ), excluding the leading particle: 114 std::map < G4int, std::vector< G4Track* > > typesAndTracks; 115 for ( std::size_t idx = 0; idx < secondariesAndPrimary.size(); ++idx ) 116 { 117 if ( idx == leadingIDX ) continue; // -- excludes the leading particle 118 auto currentTrack = secondariesAndPrimary[idx]; 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 ) 121 GROUP = -1000; // -- merge all baryons above proton/neutron in one same group -- §§ might be proposed as an option too 122 123 if ( typesAndTracks.find( GROUP ) == typesAndTracks.cend() ) 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( currentTrack ); 132 } 133 } 134 // -- and on these vectors, randomly select the surviving particles: 135 // ---- randomly select one surviving track per species 136 // ---- for this surviving track, further apply a Russian roulette 137 G4int nSecondaries = 0; // -- the number of secondaries to be returned 138 for ( auto& typeAndTrack : typesAndTracks ) 139 { 140 std::size_t nTracks = (typeAndTrack.second).size(); 141 G4Track* keptTrack; 142 // -- select one track among ones in same species: 143 if ( nTracks > 1 ) 144 { 145 auto keptTrackIDX = G4RandFlat::shootInt( nTracks ); 146 keptTrack = (typeAndTrack.second)[keptTrackIDX]; 147 keptTrack->SetWeight( keptTrack->GetWeight() * nTracks ); 148 } 149 else 150 { 151 keptTrack = (typeAndTrack.second)[0]; 152 } 153 // -- further apply a Russian Roulette on it: 154 G4bool keepTrack = false; 155 if ( fRussianRouletteKillingProbability > 0.0 ) 156 { 157 if ( G4UniformRand() > fRussianRouletteKillingProbability ) 158 { 159 keptTrack->SetWeight( keptTrack->GetWeight() / (1. - fRussianRouletteKillingProbability) ); 160 keepTrack = true; 161 } 162 } 163 else keepTrack = true; 164 if ( keepTrack ) 165 { 166 survivingMap[ keptTrack ] = true; 167 if ( keptTrack != finalStatePrimary ) ++nSecondaries; 168 } 169 } 170 // -- and if the leading is not the primary, we have to count it in nSecondaries: 171 if ( secondariesAndPrimary[leadingIDX] != finalStatePrimary ) ++nSecondaries; 172 173 // -- verify if the primary is still alive or not after above selection: 174 G4bool primarySurvived = false; 175 if ( finalStatePrimary ) 176 primarySurvived = survivingMap[ finalStatePrimary ]; 177 178 179 // -- fill the trimmed particle change: 180 // ---- fill for the primary: 181 fParticleChange.Initialize(*track); 182 if ( primarySurvived ) 183 { 184 fParticleChange.ProposeTrackStatus ( wrappedProcessParticleChange->GetTrackStatus() ); 185 fParticleChange.ProposeParentWeight ( finalStatePrimary->GetWeight() ); 186 // -- take weight from copy of primary, this one being updated in the 187 // random selection loop above 188 fParticleChange.ProposeEnergy ( finalStatePrimary->GetKineticEnergy() ); 189 fParticleChange.ProposeMomentumDirection ( finalStatePrimary->GetMomentumDirection() ); 190 } 191 else 192 { 193 fParticleChange.ProposeTrackStatus ( fStopAndKill ); 194 fParticleChange.ProposeParentWeight( 0.0 ); 195 fParticleChange.ProposeEnergy ( 0.0 ); 196 } 197 // -- fill for surviving secondaries: 198 fParticleChange.SetSecondaryWeightByProcess(true); 199 fParticleChange.SetNumberOfSecondaries(nSecondaries); 200 // ---- note we loop up to on the number of secondaries, which excludes the primary, last in secondariesAndPrimary vector: 201 ////// G4cout << callingProcess->GetProcessName() << " :"; 202 for ( auto idx = 0 ; idx < wrappedProcessParticleChange->GetNumberOfSecondaries() ; ++idx ) 203 { 204 G4Track* secondary = secondariesAndPrimary[idx]; 205 // ******************** 206 //// if ( !survivingMap[ secondary ] ) G4cout << " ["; 207 ///// else G4cout << " "; 208 ///// G4cout << secondary->GetDefinition()->GetParticleName() << " " << secondary->GetKineticEnergy(); 209 ///// if ( !survivingMap[ secondary ] ) G4cout << "]"; 210 //// if ( secondary == secondariesAndPrimary[leadingIDX] ) G4cout << " ***"; 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 change: 220 wrappedProcessParticleChange->Clear(); 221 222 if ( finalStatePrimary ) delete finalStatePrimary; 223 224 // -- finally, returns the trimmed particle change: 225 return &fParticleChange; 226 } 227