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 // particle_hp -- source file 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron transport model. 29 // 30 // 080721 To be "ClearHistories" effective, the selection scheme of angular distribution is modified 31 // by T. Koi 32 // 33 // P. Arce, Dec-2014 Conversion neutron_hp to particle_hp 34 // V. Ivanchenko, July-2023 Basic revision of particle HP classes 35 // 36 37 #include "G4ParticleHPContEnergyAngular.hh" 38 #include "G4HadronicException.hh" 39 #include "G4Neutron.hh" 40 41 G4ParticleHPContEnergyAngular::G4ParticleHPContEnergyAngular(const G4ParticleDefinition* proj) 42 { 43 theProjectile = (nullptr == proj) ? G4Neutron::Neutron() : proj; 44 currentMeanEnergy.Put(-2); 45 fCacheAngular.Put(nullptr); 46 } 47 48 G4ParticleHPContEnergyAngular::~G4ParticleHPContEnergyAngular() 49 { 50 delete[] theAngular; 51 delete fCacheAngular.Get(); 52 } 53 54 void G4ParticleHPContEnergyAngular::Init(std::istream& aDataFile) 55 { 56 aDataFile >> theTargetCode >> theAngularRep >> theInterpolation >> nEnergy; 57 const std::size_t esize = nEnergy > 0 ? nEnergy : 1; 58 theAngular = new G4ParticleHPContAngularPar[esize]; 59 theManager.Init(aDataFile); 60 for (G4int i = 0; i < nEnergy; ++i) { 61 theAngular[i].Init(aDataFile, theProjectile); 62 theAngular[i].SetInterpolation(theInterpolation); 63 #ifndef PHP_AS_HP 64 theAngular[i].PrepareTableInterpolation(); 65 #endif 66 } 67 } 68 69 G4ReactionProduct* G4ParticleHPContEnergyAngular::Sample(G4double anEnergy, G4double massCode, 70 G4double /*mass*/) 71 { 72 G4ReactionProduct* result; 73 G4int i(0); 74 G4int it(0); 75 for (i = 0; i < nEnergy; ++i) { 76 it = i; 77 #ifdef PHP_AS_HP 78 if (theAngular[i].GetEnergy() > anEnergy) break; 79 #else 80 if (theAngular[i].GetEnergy() >= anEnergy) break; 81 #endif 82 } 83 G4double targetMass = GetTarget()->GetMass(); 84 /* 85 G4cout << "G4ParticleHPContEnergyAngular::Sample E="<<anEnergy<<" code="<< massCode<< " it=" << it 86 << " M=" << targetMass << G4endl; 87 */ 88 if (it == 0) { 89 theAngular[0].SetTarget(GetTarget()); 90 theAngular[0].SetTargetCode(theTargetCode); 91 theAngular[0].SetPrimary(GetProjectileRP()); 92 result = theAngular[0].Sample(anEnergy, massCode, targetMass, theAngularRep, theInterpolation); 93 currentMeanEnergy.Put(theAngular[0].MeanEnergyOfThisInteraction()); 94 } 95 else { 96 // interpolation through alternating sampling. This needs improvement @@@ 97 // This is the cause of the He3 problem !!!!!!!! 98 // See to it, if you can improve this. 99 // G4double random = G4UniformRand(); 100 // G4double deltaE = theAngular[it].GetEnergy()-theAngular[it-1].GetEnergy(); 101 // G4double offset = theAngular[it].GetEnergy()-anEnergy; 102 // if(random<offset/deltaE) it--; 103 //--- create new 104 // if (theManager.GetScheme(0) != LINLIN) { 105 // // asserted in G4ParticleHPContEnergyAngular::init there is only one range 106 #ifdef PHP_AS_HP 107 theAngular[it].SetTarget(GetTarget()); 108 theAngular[it].SetTargetCode(theTargetCode); 109 theAngular[it].SetPrimary(GetProjectileRP()); 110 result = theAngular[it].Sample(anEnergy, massCode, targetMass, theAngularRep, theInterpolation); 111 currentMeanEnergy.Put(theAngular[it].MeanEnergyOfThisInteraction()); 112 113 #else 114 if (fCacheAngular.Get() == nullptr) { 115 auto angpar = new G4ParticleHPContAngularPar(theProjectile); 116 fCacheAngular.Put(angpar); 117 } 118 fCacheAngular.Get()->SetInterpolation(theInterpolation); 119 fCacheAngular.Get()->BuildByInterpolation(anEnergy, theManager.GetScheme(0), 120 (theAngular[it - 1]), (theAngular[it])); 121 fCacheAngular.Get()->SetTarget(GetTarget()); 122 fCacheAngular.Get()->SetTargetCode(theTargetCode); 123 fCacheAngular.Get()->SetPrimary(GetProjectileRP()); 124 125 result = 126 fCacheAngular.Get()->Sample(anEnergy, massCode, targetMass, theAngularRep, theInterpolation); 127 currentMeanEnergy.Put(fCacheAngular.Get()->MeanEnergyOfThisInteraction()); 128 #endif 129 } // end (it != 0) branch 130 //G4cout << "+!+ Emin=" << currentMeanEnergy.Get() << G4endl; 131 return result; 132 } 133 134 G4double G4ParticleHPContEnergyAngular::MeanEnergyOfThisInteraction() 135 { 136 //G4cout << "+++ Emin=" << currentMeanEnergy.Get() << G4endl; 137 G4double result(0); 138 if (currentMeanEnergy.Get() < -1.0) { 139 throw G4HadronicException(__FILE__, __LINE__, 140 "G4ParticleHPContEnergyAngular: Logical error in Product class"); 141 } 142 result = currentMeanEnergy.Get(); 143 144 currentMeanEnergy.Put(-2.0); 145 return result; 146 } 147 148 void G4ParticleHPContEnergyAngular::ClearHistories() 149 { 150 if (theAngular != nullptr) { 151 for (G4int i = 0; i < nEnergy; ++i) 152 theAngular[i].ClearHistories(); 153 } 154 155 // Added fCacheAngular ClearHistories() - this is the one actually used! 156 // Maybe theAngular does not even need ClearHistories()? 157 if (fCacheAngular.Get() != nullptr) fCacheAngular.Get()->ClearHistories(); 158 } 159