Geant4 Cross Reference |
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 // neutron_hp -- source file 26 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model. 29 // 29 // 30 // P. Arce, June-2014 Conversion neutron_hp to 30 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 31 // 31 // 32 #include "G4ParticleHPInterpolator.hh" 32 #include "G4ParticleHPInterpolator.hh" 33 33 34 #include "G4Pow.hh" 34 #include "G4Pow.hh" 35 35 36 G4double G4ParticleHPInterpolator::GetBinInteg 36 G4double G4ParticleHPInterpolator::GetBinIntegral(const G4InterpolationScheme& aScheme, 37 37 const G4double x1, const G4double x2, 38 38 const G4double y1, const G4double y2) 39 { // inline again later on @@@@ 39 { // inline again later on @@@@ 40 G4double result = 0; 40 G4double result = 0; 41 if (aScheme == HISTO || aScheme == CHISTO || 41 if (aScheme == HISTO || aScheme == CHISTO || aScheme == UHISTO) { 42 result = y1 * (x2 - x1); 42 result = y1 * (x2 - x1); 43 } 43 } 44 else if (aScheme == LINLIN || aScheme == CLI 44 else if (aScheme == LINLIN || aScheme == CLINLIN || aScheme == ULINLIN) { 45 result = 0.5 * (y2 + y1) * (x2 - x1); 45 result = 0.5 * (y2 + y1) * (x2 - x1); 46 } 46 } 47 else if (aScheme == LINLOG || aScheme == CLI 47 else if (aScheme == LINLOG || aScheme == CLINLOG || aScheme == ULINLOG) { 48 if (x1 == 0) 48 if (x1 == 0) 49 result = y1; 49 result = y1; 50 else if (x2 == 0) 50 else if (x2 == 0) 51 result = y2; 51 result = y2; 52 else { 52 else { 53 G4double b = (y2 - y1) / (G4Log(x2) - G4 53 G4double b = (y2 - y1) / (G4Log(x2) - G4Log(x1)); 54 G4double a = y1 - b * G4Log(x1); 54 G4double a = y1 - b * G4Log(x1); 55 result = (a - b) * (x2 - x1) + b * (x2 * 55 result = (a - b) * (x2 - x1) + b * (x2 * G4Log(x2) - x1 * G4Log(x1)); 56 } 56 } 57 } 57 } 58 else if (aScheme == LOGLIN || aScheme == CLO 58 else if (aScheme == LOGLIN || aScheme == CLOGLIN || aScheme == ULOGLIN) { 59 if (y1 == 0 || y2 == 0) { 59 if (y1 == 0 || y2 == 0) { 60 result = 0; 60 result = 0; 61 } 61 } 62 else { 62 else { 63 // G4double b = (std::log(y2)-std::log(y 63 // G4double b = (std::log(y2)-std::log(y1))/(x2-x1); 64 // G4double a = std::log(y1) - b*x1; 64 // G4double a = std::log(y1) - b*x1; 65 //************************************** 65 //*************************************************************** 66 // EMendoza: 66 // EMendoza: 67 // result = (std::exp(a)/b)*(std::exp(b* 67 // result = (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1)); 68 //************************************** 68 //*************************************************************** 69 if (y1 != y2) { 69 if (y1 != y2) { 70 result = (x2 - x1) * (y2 - y1) / G4Log 70 result = (x2 - x1) * (y2 - y1) / G4Log(y2 / y1); 71 } 71 } 72 else { 72 else { 73 result = y2 * (x2 - x1); 73 result = y2 * (x2 - x1); 74 } 74 } 75 //************************************** 75 //*************************************************************** 76 } 76 } 77 } 77 } 78 else if (aScheme == LOGLOG || aScheme == CLO 78 else if (aScheme == LOGLOG || aScheme == CLOGLOG || aScheme == ULOGLOG) { 79 if (x1 == 0) 79 if (x1 == 0) 80 result = y1; 80 result = y1; 81 else if (x2 == 0) 81 else if (x2 == 0) 82 result = y2; 82 result = y2; 83 else if (y1 == 0 || y2 == 0) 83 else if (y1 == 0 || y2 == 0) 84 result = 0; 84 result = 0; 85 else { 85 else { 86 G4double b = (G4Log(y2) - G4Log(y1)) / ( 86 G4double b = (G4Log(y2) - G4Log(y1)) / (G4Log(x2) - G4Log(x1)); 87 G4double a = G4Log(y1) - b * G4Log(x1); 87 G4double a = G4Log(y1) - b * G4Log(x1); 88 ; 88 ; 89 result = (G4Exp(a) / (b + 1)) 89 result = (G4Exp(a) / (b + 1)) 90 * (G4Pow::GetInstance()->powA(x 90 * (G4Pow::GetInstance()->powA(x2, b + 1) - G4Pow::GetInstance()->powA(x1, b + 1)); 91 } 91 } 92 } 92 } 93 else { 93 else { 94 throw G4HadronicException(__FILE__, __LINE 94 throw G4HadronicException(__FILE__, __LINE__, 95 "Unknown interpo 95 "Unknown interpolation scheme in G4ParticleHPVector::Integrate"); 96 } 96 } 97 return result; 97 return result; 98 } 98 } 99 G4double G4ParticleHPInterpolator::GetWeighted 99 G4double G4ParticleHPInterpolator::GetWeightedBinIntegral(const G4InterpolationScheme& aScheme, 100 100 const G4double x1, const G4double x2, 101 101 const G4double y1, const G4double y2) 102 { // inline again later on @@@@ 102 { // inline again later on @@@@ 103 G4double result = 0; 103 G4double result = 0; 104 if (aScheme == HISTO || aScheme == CHISTO || 104 if (aScheme == HISTO || aScheme == CHISTO || aScheme == UHISTO) { 105 result = 0.5 * y1 * (x2 * x2 - x1 * x1); 105 result = 0.5 * y1 * (x2 * x2 - x1 * x1); 106 } 106 } 107 else if (aScheme == LINLIN || aScheme == CLI 107 else if (aScheme == LINLIN || aScheme == CLINLIN || aScheme == ULINLIN) { 108 // G4double b = (y2-y1)/(x2-x1); 108 // G4double b = (y2-y1)/(x2-x1); 109 // G4double a = y1 - b*x1; 109 // G4double a = y1 - b*x1; 110 // result = 0.5*a*(x2*x2-x1*x1) + ( 110 // result = 0.5*a*(x2*x2-x1*x1) + (b/3.)*(x2*x2*x2-x1*x1*x1); 111 // Factor out x2-x1 to avoid divide by ze 111 // Factor out x2-x1 to avoid divide by zero 112 112 113 result = (y1 * x2 - y2 * x1) * (x2 + x1) / 113 result = (y1 * x2 - y2 * x1) * (x2 + x1) / 2. + (y2 - y1) * (x2 * x2 + x2 * x1 + x1 * x1) / 3.; 114 } 114 } 115 else if (aScheme == LINLOG || aScheme == CLI 115 else if (aScheme == LINLOG || aScheme == CLINLOG || aScheme == ULINLOG) { 116 if (x1 == 0) 116 if (x1 == 0) 117 result = y1; 117 result = y1; 118 else if (x2 == 0) 118 else if (x2 == 0) 119 result = y2; 119 result = y2; 120 else { 120 else { 121 G4double b = (y2 - y1) / (G4Log(x2) - G4 121 G4double b = (y2 - y1) / (G4Log(x2) - G4Log(x1)); 122 G4double a = y1 - b * G4Log(x1); 122 G4double a = y1 - b * G4Log(x1); 123 result = (x2 * x2 / 2. * (a - b / 2. + b 123 result = (x2 * x2 / 2. * (a - b / 2. + b * G4Log(x2))) 124 - (x1 * x1 / 2. * (a - b / 2. + 124 - (x1 * x1 / 2. * (a - b / 2. + b * G4Log(x1))); 125 } 125 } 126 } 126 } 127 else if (aScheme == LOGLIN || aScheme == CLO 127 else if (aScheme == LOGLIN || aScheme == CLOGLIN || aScheme == ULOGLIN) { 128 if (y1 == 0 || y2 == 0) 128 if (y1 == 0 || y2 == 0) 129 result = 0; 129 result = 0; 130 else { 130 else { 131 G4double b = (G4Log(y2) - G4Log(y1)) / ( 131 G4double b = (G4Log(y2) - G4Log(y1)) / (x2 - x1); 132 G4double a = G4Log(y1) - b * x1; 132 G4double a = G4Log(y1) - b * x1; 133 result = G4Exp(a) / (b * b) * (G4Exp(b * 133 result = G4Exp(a) / (b * b) * (G4Exp(b * x2) * (b * x2 - 1.) - G4Exp(b * x1) * (b * x1 - 1.)); 134 } 134 } 135 } 135 } 136 else if (aScheme == LOGLOG || aScheme == CLO 136 else if (aScheme == LOGLOG || aScheme == CLOGLOG || aScheme == ULOGLOG) { 137 if (x1 == 0) 137 if (x1 == 0) 138 result = y1; 138 result = y1; 139 else if (x2 == 0) 139 else if (x2 == 0) 140 result = y2; 140 result = y2; 141 else if (y1 == 0 || y2 == 0) 141 else if (y1 == 0 || y2 == 0) 142 result = 0; 142 result = 0; 143 else { 143 else { 144 G4double b = (G4Log(y2) - G4Log(y1)) / ( 144 G4double b = (G4Log(y2) - G4Log(y1)) / (G4Log(x2) - G4Log(x1)); 145 G4double a = G4Log(y1) - b * G4Log(x1); 145 G4double a = G4Log(y1) - b * G4Log(x1); 146 ; 146 ; 147 result = G4Exp(a) / (b + 2.) 147 result = G4Exp(a) / (b + 2.) 148 * (G4Pow::GetInstance()->powA(x 148 * (G4Pow::GetInstance()->powA(x2, b + 2.) - G4Pow::GetInstance()->powA(x1, b + 2)); 149 } 149 } 150 } 150 } 151 else { 151 else { 152 throw G4HadronicException(__FILE__, __LINE 152 throw G4HadronicException(__FILE__, __LINE__, 153 "Unknown interpo 153 "Unknown interpolation scheme in G4ParticleHPVector::Integrate"); 154 } 154 } 155 return result; 155 return result; 156 } 156 } 157 157