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 ]

Diff markup

Differences between /processes/hadronic/models/lend/src/G4LENDInelastic.cc (Version 11.3.0) and /processes/hadronic/models/lend/src/G4LENDInelastic.cc (Version 9.5)


  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 ////////////////////////////////////////////// << 
 27 //                                             << 
 28 //  File:   G4LENDInelastic.cc                 << 
 29 //  Date:   24 March 2020                      << 
 30 //  Author: Dennis Wright                      << 
 31 //                                             << 
 32 //  Description: model for inelastic scatterin << 
 33 //               gammas at energies of order 2 << 
 34 //               This model uses GIDI particle << 
 35 //               in spectrum mode.  In this mo << 
 36 //               for each possible particle ty << 
 37 //               given interaction.  Unlike Ge << 
 38 //               consideration of event-by-eve << 
 39 //               Indeed, forcing such conserva << 
 40 //               introduces correlations and d << 
 41 //               spectra which are not present << 
 42 //                                             << 
 43 //               In order to use GIDI data wit << 
 44 //               minimal event-by-event baryon << 
 45 //               enforced which allows deviati << 
 46 //               giving warnings.  Neither cha << 
 47 //               conservation is enforced.  Un << 
 48 //               fragment (n, p, d, t, alpha)  << 
 49 //               after a large number of event << 
 50 //               conservation also approach th << 
 51 //               this limit.  The mass, charge << 
 52 //               large fragments, however, are << 
 53 //               data very well.  This is a re << 
 54 //               baryon number conservation an << 
 55 //               fragment spectra are correct. << 
 56 //                                             << 
 57 ////////////////////////////////////////////// << 
 58                                                << 
 59 #include "G4LENDInelastic.hh"                      26 #include "G4LENDInelastic.hh"
 60 #include "G4SystemOfUnits.hh"                  <<  27 
 61 #include "G4Nucleus.hh"                            28 #include "G4Nucleus.hh"
 62 #include "G4IonTable.hh"                       <<  29 #include "G4ParticleTable.hh"
 63 #include <algorithm>                           << 
 64 #include <random>                              << 
 65                                                    30   
 66 G4HadFinalState* G4LENDInelastic::ApplyYoursel <<  31 G4HadFinalState * G4LENDInelastic::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& aTarg )
 67                                                << 
 68 {                                                  32 {
 69   G4ThreeVector projMom = aTrack.Get4Momentum( << 
 70   G4double temp = aTrack.GetMaterial()->GetTem << 
 71                                                    33 
 72   G4int iZ = aTarg.GetZ_asInt();               <<  34    G4ThreeVector proj_p = aTrack.Get4Momentum().vect();
 73   G4int iA = aTarg.GetA_asInt();               <<  35 
 74   G4int iM = 0;                                <<  36    G4double temp = aTrack.GetMaterial()->GetTemperature();
 75   if (aTarg.GetIsotope() != nullptr) iM = aTar <<  37 
 76                                                <<  38    //G4int iZ = int ( aTarg.GetZ() );
 77   G4double ke = aTrack.GetKineticEnergy();     <<  39    //G4int iA = int ( aTarg.GetN() );
 78                                                <<  40    //migrate to integer A and Z (GetN_asInt returns number of neutrons in the nucleus since this) 
 79   G4HadFinalState* theResult = &theParticleCha <<  41    G4int iZ = aTarg.GetZ_asInt();
 80   theResult->Clear();                          <<  42    G4int iA = aTarg.GetA_asInt();
 81                                                <<  43    //G4cout << "target: Z = " << iZ << " N = " << iA << G4endl;
 82   G4GIDI_target* aGIDITarget =                 <<  44 
 83     get_target_from_map(lend_manager->GetNucle <<  45    G4double ke = aTrack.GetKineticEnergy();
 84   if (aGIDITarget == nullptr) {                <<  46    //G4cout << "projectile: KE = " << ke/MeV << " [MeV]" << G4endl;
 85 //    G4cout << " No target found " << G4endl; <<  47 
 86     theParticleChange.Clear();                 <<  48    G4HadFinalState* theResult = &theParticleChange;
 87     theParticleChange.SetStatusChange(isAlive) <<  49    theResult->Clear();
 88     theParticleChange.SetEnergyChange(aTrack.G <<  50 
 89     theParticleChange.SetMomentumChange(aTrack <<  51    G4GIDI_target* aTarget = usedTarget_map.find( lend_manager->GetNucleusEncoding( iZ , iA ) )->second->GetTarget();
 90     return &theParticleChange;                 <<  52    std::vector<G4GIDI_Product>* products = aTarget->getOthersFinalState( ke*MeV, temp, NULL, NULL );
 91   }                                            <<  53    if ( products != NULL ) 
 92 // return returnUnchanged(aTrack, theResult);  <<  54    {
 93                                                <<  55 
 94   // Get GIDI final state products for givern  <<  56       G4ThreeVector psum(0);
 95   G4int loop(0);                               <<  57 
 96   G4int loopMax = 1000;                        <<  58       int totN = 0;
 97   std::vector<G4GIDI_Product>* products;       <<  59       for ( G4int j = 0; j < int( products->size() ); j++ ) 
 98   do {                                         <<  60       {
 99     products = aGIDITarget->getOthersFinalStat <<  61 
100     loop++;                                    <<  62          G4int jZ = (*products)[j].Z; 
101   } while (products == nullptr && loop < loopM <<  63          G4int jA = (*products)[j].A; 
102                                                <<  64 
103   // G4LENDInelastic accepts all light fragmen <<  65          //G4cout << "ZA = " << 1000 * (*products)[j].Z + (*products)[j].A << "  EK = "
104   // and removes any heavy fragments which cau <<  66          //     << (*products)[j].kineticEnergy
105   // Charge and energy non-conservation still  <<  67          //     << " px  " <<  (*products)[j].px
106   // of events, this improves on average.      <<  68          //     << " py  " <<  (*products)[j].py
107                                                <<  69          //     << " pz  " <<  (*products)[j].pz
108   if (loop > loopMax - 1) {                    <<  70          //     << G4endl;
109 //    G4cout << " too many loops, return initi <<  71 
110                                                <<  72          G4DynamicParticle* theSec = new G4DynamicParticle;
111     theParticleChange.Clear();                 <<  73 
112     theParticleChange.SetStatusChange(isAlive) <<  74          if ( jA == 1 && jZ == 1 )
113     theParticleChange.SetEnergyChange(aTrack.G <<  75          {
114     theParticleChange.SetMomentumChange(aTrack <<  76             theSec->SetDefinition( G4Proton::Proton() );
115     return &theParticleChange;                 <<  77             totN += 1;
116                                                <<  78          }
117 //    if (aTrack.GetDefinition() == G4Proton:: <<  79          else if ( jA == 1 && jZ == 0 )
118 //        aTrack.GetDefinition() == G4Neutron: <<  80          {
119 //       theResult = preco->ApplyYourself(aTra <<  81             theSec->SetDefinition( G4Neutron::Neutron() );
120 //     } else {                                <<  82             totN += 1;
121 //       theResult = returnUnchanged(aTrack, t <<  83          } 
122 //     }                                       <<  84          else if ( jZ > 0 )
123                                                <<  85          {
124   } else {                                     <<  86             if ( jA != 0 )
125     G4int iTotZ = iZ + aTrack.GetDefinition()- <<  87             {
126     G4int iTotA = iA + aTrack.GetDefinition()- <<  88                theSec->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( jZ , jA , 0 , 0 ) );
127                                                <<  89                totN += jA;
128     // Loop over GIDI products and separate li <<  90             }
129     G4int GZtot(0);                            <<  91             else 
130     G4int GAtot(0);                            <<  92             {
131     G4int productA(0);                         <<  93                theSec->SetDefinition( G4ParticleTable::GetParticleTable()->FindIon( jZ , iA+1-totN , 0 , 0 ) );
132     G4int productZ(0);                         <<  94             }
133     std::vector<G4int> lightProductIndex;      <<  95          } 
134     std::vector<G4int> heavyProductIndex;      <<  96          else
135     for (G4int i = 0; i < int( products->size( <<  97          {
136       productA = (*products)[i].A;             <<  98             theSec->SetDefinition( G4Gamma::Gamma() );
137       if (productA < 5) {                      <<  99          } 
138         lightProductIndex.push_back(i);        << 100 
139         GZtot += (*products)[i].Z;             << 101          G4ThreeVector p( (*products)[j].px*MeV , (*products)[j].py*MeV , (*products)[j].pz*MeV ); 
140         GAtot += productA;                     << 102          psum += p; 
141       } else {                                 << 103          if ( p.mag() == 0 ) p = proj_p - psum;
142         heavyProductIndex.push_back(i);        << 104 
143       }                                        << 105          theSec->SetMomentum( p );
144     }                                          << 106 
145                                                << 107          theResult->AddSecondary( theSec );
146     // Randomize order of heavies to correct s << 108       } 
147     // std::random_shuffle(heavyProductIndex.b << 109    }
148     // std::cout << " Heavy product index befo << 110    delete products;
149     // for (G4int i = 0; i < int(heavyProductI << 111 
150     // std::cout << std::endl;                 << 112    theResult->SetStatusChange( stopAndKill );
151                                                << 113 
152     auto rng = std::default_random_engine {};  << 114    return theResult; 
153     std::shuffle(heavyProductIndex.begin(), he << 115 
154                                                << 
155     // std::cout << " Heavy product index afte << 
156     // for (G4int i = 0; i < int(heavyProductI << 
157     // std::cout << std::endl;                 << 
158                                                << 
159     std::vector<G4int> savedHeavyIndex;        << 
160     G4int itest(0);                            << 
161     for (G4int i = 0; i < int(heavyProductInde << 
162       itest = heavyProductIndex[i];            << 
163       productA = (*products)[itest].A;         << 
164       productZ = (*products)[itest].Z;         << 
165       if ((GAtot + productA <= iTotA) && (GZto << 
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(lightProductInde << 
175       itest = lightProductIndex[k];            << 
176       G4cout << "(" << (*products)[itest].Z << << 
177     }                                          << 
178     G4cout << G4endl;                          << 
179                                                << 
180     G4cout << " saved heavy products = ";      << 
181     for (G4int k = 0; k < int(savedHeavyIndex. << 
182       itest = savedHeavyIndex[k];              << 
183       G4cout << "(" << (*products)[itest].Z << << 
184     }                                          << 
185     G4cout << G4endl;                          << 
186 */                                             << 
187     // Now convert saved products to Geant4 pa << 
188     // Note that, at least for heavy fragments << 
189     // have slightly different values.         << 
190                                                << 
191     G4DynamicParticle* theSec = nullptr;       << 
192     G4ThreeVector Psum;                        << 
193     for (G4int i = 0; i < int(lightProductInde << 
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::Neutr << 
200       } else if (productA == 1 && productZ ==  << 
201         theSec->SetDefinition(G4Proton::Proton << 
202       } else if (productA == 2 && productZ ==  << 
203         theSec->SetDefinition(G4Deuteron::Deut << 
204       } else if (productA == 3 && productZ ==  << 
205         theSec->SetDefinition(G4Triton::Triton << 
206       } else if (productA == 4 && productZ ==  << 
207         theSec->SetDefinition(G4Alpha::Alpha() << 
208       } else {                                 << 
209         theSec->SetDefinition(G4Gamma::Gamma() << 
210       }                                        << 
211                                                << 
212       G4ThreeVector momentum((*products)[itest << 
213                              (*products)[itest << 
214                              (*products)[itest << 
215       Psum += momentum;                        << 
216       theSec->SetMomentum(momentum);           << 
217 //      theResult->AddSecondary(theSec);       << 
218       theParticleChange.AddSecondary(theSec, s << 
219     }                                          << 
220                                                << 
221     G4int productM(0);                         << 
222     for (G4int i = 0; i < int(savedHeavyIndex. << 
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::GetIon << 
229                                                << 
230                                                << 
231       G4ThreeVector momentum((*products)[itest << 
232                              (*products)[itest << 
233                              (*products)[itest << 
234       Psum += momentum;                        << 
235       theSec->SetMomentum(momentum);           << 
236 //      theResult->AddSecondary(theSec);       << 
237       theParticleChange.AddSecondary(theSec, s << 
238     }                                          << 
239                                                << 
240     // Create heavy fragment if necessary to t << 
241     // Note: this step is only required to pre << 
242     //       where "catastrophic" non-conserva << 
243     //       The residual generated will not n << 
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 << 
249         // Violate charge conservation and set << 
250         // G4cout << " Z = 1, A = "<< iTotA -  << 
251         theSec->SetDefinition(G4IonTable::GetI << 
252       } else {                                 << 
253         theSec->SetDefinition(G4IonTable::GetI << 
254       }                                        << 
255       theSec->SetMomentum(projMom - Psum);     << 
256 //      theResult->AddSecondary(theSec);       << 
257       theParticleChange.AddSecondary(theSec, s << 
258     }                                          << 
259   } // loop OK                                 << 
260                                                << 
261   delete products;                             << 
262 //  theResult->SetStatusChange( stopAndKill ); << 
263   theParticleChange.SetStatusChange( stopAndKi << 
264 //  return theResult;                          << 
265   return &theParticleChange;                   << 
266 }                                                 116 }
267                                                   117