ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LENDInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LENDInelastic.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 //
27 // //
28 // File: G4LENDInelastic.cc //
29 // Date: 24 March 2020 //
30 // Author: Dennis Wright //
31 // //
32 // Description: model for inelastic scattering of neutrons, light ions and //
33 // gammas at energies of order 20 MeV and lower. //
34 // This model uses GIDI particle data which are stored mostly //
35 // in spectrum mode. In this mode, spectra are reproduced //
36 // for each possible particle type which can result from a //
37 // given interaction. Unlike Geant4, this is done without //
38 // consideration of event-by-event conservation rules. //
39 // Indeed, forcing such conservation on GIDI output products //
40 // introduces correlations and distortions in the resulting //
41 // spectra which are not present in the data. //
42 // //
43 // In order to use GIDI data within the Geant4 framework, a //
44 // minimal event-by-event baryon number conservation is //
45 // enforced which allows deviations of up to 1 GeV without //
46 // giving warnings. Neither charge, nor energy, nor momentum //
47 // conservation is enforced. Under this scheme, light //
48 // fragment (n, p, d, t, alpha) spectra are well reproduced //
49 // after a large number of events. Charge and energy //
50 // conservation also approach their event-by-event values in //
51 // this limit. The mass, charge and energy distributions of //
52 // large fragments, however, are not expected to reproduce the //
53 // data very well. This is a result of forcing the crude //
54 // baryon number conservation and ensuring that the light //
55 // fragment spectra are correct. //
56 // //
58 
59 #include "G4LENDInelastic.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "G4Nucleus.hh"
62 #include "G4IonTable.hh"
63 #include <algorithm>
64 #include <random>
65 
67  G4Nucleus& aTarg)
68 {
69  G4ThreeVector projMom = aTrack.Get4Momentum().vect();
70  G4double temp = aTrack.GetMaterial()->GetTemperature();
71 
72  G4int iZ = aTarg.GetZ_asInt();
73  G4int iA = aTarg.GetA_asInt();
74  G4int iM = 0;
75  if (aTarg.GetIsotope() != nullptr) iM = aTarg.GetIsotope()->Getm();
76 
77  G4double ke = aTrack.GetKineticEnergy();
78 
79  G4HadFinalState* theResult = &theParticleChange;
80  theResult->Clear();
81 
82  G4GIDI_target* aGIDITarget =
84  if (aGIDITarget == nullptr) {
85 // G4cout << " No target found " << G4endl;
90  return &theParticleChange;
91  }
92 // return returnUnchanged(aTrack, theResult);
93 
94  // Get GIDI final state products for givern target and projectile
95  G4int loop(0);
96  G4int loopMax = 1000;
97  std::vector<G4GIDI_Product>* products;
98  do {
99  products = aGIDITarget->getOthersFinalState(ke*MeV, temp, MyRNG, NULL);
100  loop++;
101  } while (products == nullptr && loop < loopMax);
102 
103  // G4LENDInelastic accepts all light fragments and gammas from GIDI (A < 5)
104  // and removes any heavy fragments which cause large baryon number violation.
105  // Charge and energy non-conservation still occur, but over a large number
106  // of events, this improves on average.
107 
108  if (loop > loopMax - 1) {
109 // G4cout << " too many loops, return intial state " << G4endl;
110 
115  return &theParticleChange;
116 
117 // if (aTrack.GetDefinition() == G4Proton::Proton() ||
118 // aTrack.GetDefinition() == G4Neutron::Neutron() ) {
119 // theResult = preco->ApplyYourself(aTrack, aTarg);
120 // } else {
121 // theResult = returnUnchanged(aTrack, theResult);
122 // }
123 
124  } else {
125  G4int iTotZ = iZ + aTrack.GetDefinition()->GetAtomicNumber();
126  G4int iTotA = iA + aTrack.GetDefinition()->GetAtomicMass();
127 
128  // Loop over GIDI products and separate light from heavy fragments
129  G4int GZtot(0);
130  G4int GAtot(0);
131  G4int productA(0);
132  G4int productZ(0);
133  std::vector<G4int> lightProductIndex;
134  std::vector<G4int> heavyProductIndex;
135  for (G4int i = 0; i < int( products->size() ); i++ ) {
136  productA = (*products)[i].A;
137  if (productA < 5) {
138  lightProductIndex.push_back(i);
139  GZtot += (*products)[i].Z;
140  GAtot += productA;
141  } else {
142  heavyProductIndex.push_back(i);
143  }
144  }
145 
146  // Randomize order of heavies to correct somewhat for sampling bias
147  // std::random_shuffle(heavyProductIndex.begin(), heavyProductIndex.end() );
148  // std::cout << " Heavy product index before shuffle : " ;
149  // for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) std::cout << heavyProductIndex[i] << ", " ;
150  // std::cout << std::endl;
151 
152  auto rng = std::default_random_engine {};
153  std::shuffle(heavyProductIndex.begin(), heavyProductIndex.end(), rng);
154 
155  // std::cout << " Heavy product index after shuffle : " ;
156  // for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) std::cout << heavyProductIndex[i] << ", " ;
157  // std::cout << std::endl;
158 
159  std::vector<G4int> savedHeavyIndex;
160  G4int itest(0);
161  for (G4int i = 0; i < int(heavyProductIndex.size() ); i++) {
162  itest = heavyProductIndex[i];
163  productA = (*products)[itest].A;
164  productZ = (*products)[itest].Z;
165  if ((GAtot + productA <= iTotA) && (GZtot + productZ <= iTotZ) ) {
166  savedHeavyIndex.push_back(itest);
167  GZtot += productZ;
168  GAtot += productA;
169  }
170  }
171 
172 /*
173  G4cout << " saved light products = ";
174  for (G4int k = 0; k < int(lightProductIndex.size() ); k++ ) {
175  itest = lightProductIndex[k];
176  G4cout << "(" << (*products)[itest].Z << ", " << (*products)[itest].A << "), ";
177  }
178  G4cout << G4endl;
179 
180  G4cout << " saved heavy products = ";
181  for (G4int k = 0; k < int(savedHeavyIndex.size() ); k++ ) {
182  itest = savedHeavyIndex[k];
183  G4cout << "(" << (*products)[itest].Z << ", " << (*products)[itest].A << "), ";
184  }
185  G4cout << G4endl;
186 */
187  // Now convert saved products to Geant4 particles
188  // Note that, at least for heavy fragments, GIDI masses and Geant4 masses
189  // have slightly different values.
190 
191  G4DynamicParticle* theSec = nullptr;
192  G4ThreeVector Psum;
193  for (G4int i = 0; i < int(lightProductIndex.size() ); i++) {
194  itest = lightProductIndex[i];
195  productZ = (*products)[itest].Z;
196  productA = (*products)[itest].A;
197  theSec = new G4DynamicParticle();
198  if (productA == 1 && productZ == 0) {
199  theSec->SetDefinition(G4Neutron::Neutron() );
200  } else if (productA == 1 && productZ == 1) {
201  theSec->SetDefinition(G4Proton::Proton() );
202  } else if (productA == 2 && productZ == 1) {
204  } else if (productA == 3 && productZ == 1) {
205  theSec->SetDefinition(G4Triton::Triton() );
206  } else if (productA == 4 && productZ == 2) {
207  theSec->SetDefinition(G4Alpha::Alpha() );
208  } else {
209  theSec->SetDefinition(G4Gamma::Gamma() );
210  }
211 
212  G4ThreeVector momentum((*products)[itest].px*MeV,
213  (*products)[itest].py*MeV,
214  (*products)[itest].pz*MeV );
215  Psum += momentum;
216  theSec->SetMomentum(momentum);
217 // theResult->AddSecondary(theSec);
219  }
220 
221  G4int productM(0);
222  for (G4int i = 0; i < int(savedHeavyIndex.size() ); i++) {
223  itest = savedHeavyIndex[i];
224  productZ = (*products)[itest].Z;
225  productA = (*products)[itest].A;
226  productM = (*products)[itest].m;
227  theSec = new G4DynamicParticle();
228  theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(productZ,
229  productA,
230  productM) );
231  G4ThreeVector momentum((*products)[itest].px*MeV,
232  (*products)[itest].py*MeV,
233  (*products)[itest].pz*MeV );
234  Psum += momentum;
235  theSec->SetMomentum(momentum);
236 // theResult->AddSecondary(theSec);
238  }
239 
240  // Create heavy fragment if necessary to try to balance A, Z
241  // Note: this step is only required to prevent warnings at the process level
242  // where "catastrophic" non-conservation tolerances are set to 1 GeV.
243  // The residual generated will not necessarily be the one that would
244  // occur in the actual reaction.
245  if (iTotA - GAtot > 1) {
246  theSec = new G4DynamicParticle();
247  if (iTotZ == GZtot) {
248  // Special case when a nucleus of only neutrons is requested
249  // Violate charge conservation and set Z = 1
250  // G4cout << " Z = 1, A = "<< iTotA - GAtot << " created " << G4endl;
251  theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(1, iTotA-GAtot, 0) );
252  } else {
253  theSec->SetDefinition(G4IonTable::GetIonTable()->GetIon(iTotZ-GZtot, iTotA-GAtot, 0) );
254  }
255  theSec->SetMomentum(projMom - Psum);
256 // theResult->AddSecondary(theSec);
258  }
259  } // loop OK
260 
261  delete products;
262 // theResult->SetStatusChange( stopAndKill );
264 // return theResult;
265  return &theParticleChange;
266 }