Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/G4LENDInelastic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 ////////////////////////////////////////////////////////////////////////////////
 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 //                                                                            //
 57 ////////////////////////////////////////////////////////////////////////////////
 58 
 59 #include "G4LENDInelastic.hh"
 60 #include "G4SystemOfUnits.hh"
 61 #include "G4Nucleus.hh"
 62 #include "G4IonTable.hh"
 63 #include <algorithm>
 64 #include <random>
 65   
 66 G4HadFinalState* G4LENDInelastic::ApplyYourself(const G4HadProjectile& aTrack,
 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 =
 83     get_target_from_map(lend_manager->GetNucleusEncoding(iZ, iA, iM) );
 84   if (aGIDITarget == nullptr) {
 85 //    G4cout << " No target found " << G4endl; 
 86     theParticleChange.Clear();
 87     theParticleChange.SetStatusChange(isAlive);
 88     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
 89     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
 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 initial state " << G4endl;
110 
111     theParticleChange.Clear();
112     theParticleChange.SetStatusChange(isAlive);
113     theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
114     theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
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) {
203         theSec->SetDefinition(G4Deuteron::Deuteron() );
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);
218       theParticleChange.AddSecondary(theSec, secID);
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);
237       theParticleChange.AddSecondary(theSec, secID);
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);
257       theParticleChange.AddSecondary(theSec, secID);
258     }
259   } // loop OK
260 
261   delete products;
262 //  theResult->SetStatusChange( stopAndKill );
263   theParticleChange.SetStatusChange( stopAndKill );
264 //  return theResult; 
265   return &theParticleChange;
266 }
267