ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BOptnLeadingParticle.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4BOptnLeadingParticle.cc
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 //
28 
29 #include <vector>
30 #include <map>
31 
32 
34  : G4VBiasingOperation ( name ),
35  fRussianRouletteKillingProbability ( -1.0 )
36 {
37 }
38 
40 {
41 }
42 
44  const G4Track* track,
45  const G4Step* step,
46  G4bool& )
47 {
48  // -- collect wrapped process particle change:
49  auto wrappedProcessParticleChange = callingProcess->GetWrappedProcess()->PostStepDoIt(*track,*step);
50 
51  // -- does nothing in case the primary stays alone or in weird situation where all are killed...
52  if ( wrappedProcessParticleChange->GetNumberOfSecondaries() == 0 ) return wrappedProcessParticleChange;
53  if ( wrappedProcessParticleChange->GetTrackStatus() == fKillTrackAndSecondaries ) return wrappedProcessParticleChange;
54  // -- ... else, collect the secondaries in a same vector (the primary is pushed in this vector, if surviving, later see [**]):
55  std::vector < G4Track* > secondariesAndPrimary;
56  for ( auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary.push_back( wrappedProcessParticleChange->GetSecondary( i ) );
57 
58 
59  // -- If case the primary survives, need to collect its new state. In the general case of the base class G4VParticleChange
60  // -- this is tricky, as this class does not hold the primary changes (and we have to build a fake step and fake track
61  // -- caring about the touchables, etc.) So we limit here to the G4ParticleChange case, checking the reality of this
62  // -- class with a dynamic cast. If we don't have such an actual G4DynamicParticle object, we give up the biasing and return
63  // -- the untrimmed process final state.
64  // -- Note that this case does not happen often, as the technique is intended for inelastic processes. For case where several
65  // -- particles can be produced without killing the primary, we have for example the electron-/positron-nuclear
66  G4ParticleChange* castParticleChange ( nullptr );
67  G4Track* finalStatePrimary ( nullptr );
68  if ( ( wrappedProcessParticleChange->GetTrackStatus() != fStopAndKill ) )
69  {
70  // fFakePrimaryTrack->CopyTrackInfo( *track );
71  // fFakeStep ->InitializeStep( fFakePrimaryTrack );
72  // wrappedProcessParticleChange->UpdateStepForPostStep( fFakeStep );
73  // fFakeStep->UpdateTrack();
74  castParticleChange = dynamic_cast< G4ParticleChange* >( wrappedProcessParticleChange );
75  if ( castParticleChange == nullptr )
76  {
77  G4cout << " **** G4BOptnLeadingParticle::ApplyFinalStateBiasing(...) : can not bias for " << callingProcess->GetProcessName() << ", this is just a warning." << G4endl;
78  return wrappedProcessParticleChange;
79  }
80  finalStatePrimary = new G4Track( *track );
81  finalStatePrimary->SetKineticEnergy ( castParticleChange->GetEnergy() );
82  finalStatePrimary->SetWeight ( castParticleChange->GetWeight() );
83  finalStatePrimary->SetMomentumDirection( *(castParticleChange->GetMomentumDirection()) );
84  // -- [**] push the primary as the last track in the vector of tracks:
85  secondariesAndPrimary.push_back( finalStatePrimary );
86  }
87 
88  // -- Ensure the secondaries all have the primary weight:
89  // ---- collect primary track weight, from updated by process if alive, or from original copy if died:
90  G4double primaryWeight;
91  if ( finalStatePrimary ) primaryWeight = finalStatePrimary->GetWeight();
92  else primaryWeight = track ->GetWeight();
93  // ---- now set this same weight to all secondaries:
94  for ( auto i = 0 ; i < wrappedProcessParticleChange->GetNumberOfSecondaries() ; i++ ) secondariesAndPrimary[ i ]->SetWeight( primaryWeight );
95 
96 
97  // -- finds the leading particle, initialize a map of surviving tracks, tag as surviving the leading track:
98  size_t leadingIDX = 0;
99  G4double leadingEnergy = -1;
100  std::map< G4Track*, G4bool > survivingMap;
101  for ( size_t idx = 0; idx < secondariesAndPrimary.size(); idx++ )
102  {
103  survivingMap[ secondariesAndPrimary[idx] ] = false;
104  if ( secondariesAndPrimary[idx]->GetKineticEnergy() > leadingEnergy )
105  {
106  leadingEnergy = secondariesAndPrimary[idx]->GetKineticEnergy();
107  leadingIDX = idx;
108  }
109  }
110  survivingMap[ secondariesAndPrimary[leadingIDX] ] = true; // -- tag as surviving the leading particle
111 
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 ( 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 ) 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  }
132  }
133  // -- and on these vectors, randomly select the surviving particles:
134  // ---- randomly select one surviving track per species
135  // ---- for this surviving track, further apply a Russian roulette
136  G4int nSecondaries = 0; // -- the number of secondaries to be returned
137  for ( auto typeAndTrack : typesAndTracks )
138  {
139  size_t nTracks = (typeAndTrack.second).size();
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;
155  {
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  }
169  // -- and if the leading is not the primary, we have to count it in nSecondaries:
170  if ( secondariesAndPrimary[leadingIDX] != finalStatePrimary ) nSecondaries++;
171 
172  // -- verify if the primary is still alive or not after above selection:
173  G4bool primarySurvived = false;
174  if ( finalStatePrimary ) primarySurvived = survivingMap[ finalStatePrimary ];
175 
176 
177  // -- fill the trimmed particle change:
178  // ---- fill for the primary:
179  fParticleChange.Initialize(*track);
180  if ( primarySurvived )
181  {
182  fParticleChange.ProposeTrackStatus ( wrappedProcessParticleChange->GetTrackStatus() );
183  fParticleChange.ProposeParentWeight ( finalStatePrimary->GetWeight() ); // -- take weight from copy of primary, this one being updated in the random selection loop above
184  fParticleChange.ProposeEnergy ( finalStatePrimary->GetKineticEnergy() );
186  }
187  else
188  {
192  }
193  // -- fill for surviving secondaries:
196  // ---- note we loop up to on the number of secondaries, which excludes the primary, last in secondariesAndPrimary vector:
198  for ( auto idx = 0 ; idx < wrappedProcessParticleChange->GetNumberOfSecondaries() ; idx++ )
199  {
200  G4Track* secondary = secondariesAndPrimary[idx];
201  // ********************
207  // ******************
208  if ( survivingMap[ secondary ] ) fParticleChange.AddSecondary( secondary );
209  else delete secondary;
210  }
212 
213  // -- clean the wrapped process particle change:
214  wrappedProcessParticleChange->Clear();
215 
216  if ( finalStatePrimary ) delete finalStatePrimary;
217 
218  // -- finally, returns the trimmed particle change:
219  return &fParticleChange;
220 
221 }