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 // $Id: G4SauterGavrilaAngularDistribution.cc 68797 2013-04-05 13:27:11Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class file 30 // GEANT4 Class file 30 // 31 // 31 // 32 // 32 // File name: G4SauterGavrilaAngularDistri 33 // File name: G4SauterGavrilaAngularDistribution 33 // 34 // 34 // Author: Vladimir Ivanchenko using Miche 35 // Author: Vladimir Ivanchenko using Michel Maire algorithm 35 // developed for Geant3 36 // developed for Geant3 36 // 37 // 37 // Creation date: 23 July 2012 38 // Creation date: 23 July 2012 38 // 39 // 39 // Modified: << 40 // 40 // 04.05.2017 Marilena Bandieramonte impleme << 41 // ------------------------------------------- 41 // ------------------------------------------------------------------- 42 // 42 // 43 43 44 #include "G4SauterGavrilaAngularDistribution.h 44 #include "G4SauterGavrilaAngularDistribution.hh" 45 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh" 46 #include "Randomize.hh" 46 #include "Randomize.hh" 47 47 48 G4SauterGavrilaAngularDistribution::G4SauterGa 48 G4SauterGavrilaAngularDistribution::G4SauterGavrilaAngularDistribution() 49 : G4VEmAngularDistribution("SauterGavrila") << 49 : G4VEmAngularDistribution("AngularGenSauterGavrila") 50 {} 50 {} 51 51 52 G4SauterGavrilaAngularDistribution::~G4SauterG << 52 G4SauterGavrilaAngularDistribution::~G4SauterGavrilaAngularDistribution() 53 = default; << 53 {} 54 54 55 G4ThreeVector& G4SauterGavrilaAngularDistribu << 55 G4ThreeVector& >> 56 G4SauterGavrilaAngularDistribution::SampleDirection( 56 const G4DynamicParticle* dp, G4double, 57 const G4DynamicParticle* dp, G4double, G4int, const G4Material*) 57 { 58 { 58 static const G4double emin = 1*CLHEP::eV; << 59 G4double tau = dp->GetKineticEnergy()/electron_mass_c2; 59 static const G4double emax = 100*CLHEP::MeV; << 60 static const G4double taulimit = 30.0; 60 61 61 G4double energy = std::max(dp->GetKineticEne << 62 if (tau > taulimit) { 62 if (energy > emax) { << 63 fLocalDirection = dp->GetMomentumDirection(); 63 fLocalDirection = dp->GetMomentumDirection << 64 // Bugzilla 1120 >> 65 // SI on 05/09/2010 as suggested by JG 04/09/10 64 } else { 66 } else { 65 // Initial algorithm according Penelope 20 << 67 66 // F.Sauter Ann. Physik 9, 217(1931); 11, << 68 G4double invgamma = 1.0/(tau + 1.0); 67 // Modified according Penelope 2014 manual << 69 G4double beta = std::sqrt(tau*(tau + 2.0))*invgamma; 68 G4double costheta = 0.0; << 70 G4double b = 0.5*tau*(tau*tau - 1.0); 69 << 71 G4double invgamma2 = invgamma*invgamma; 70 // 1) initialize energy-dependent variable << 72 71 // Variable naming according to Eq. (2.24) << 73 G4double rndm,term,greject,grejsup,costeta,sint2; 72 // (pag. 44) << 74 if (tau < 1.) { grejsup = (1.+b-beta*b)/invgamma2; } 73 G4double tau = energy/electron_mass_c2; << 75 else { grejsup = (1.+b+beta*b)/invgamma2; } 74 G4double gamma = 1.0 + tau; << 76 75 G4double beta = std::sqrt(tau*(tau + 2.0)) << 77 do { 76 << 78 rndm = 1 - 2*G4UniformRand(); 77 // ac corresponds to "A" of Eq. (2.31) << 79 costeta = (rndm + beta)/(rndm*beta + 1); 78 // << 80 term = invgamma2/(1 + beta*rndm); 79 G4double ac = (1.0 - beta)/beta; << 81 sint2 = (1 - costeta)*(1 + costeta); 80 G4double a1 = 0.5*beta*gamma*tau*(gamma-2. << 82 greject = sint2*(1 + b*term)/(term*term); 81 G4double a2 = ac + 2.0; << 83 82 // gtmax = maximum of the rejection functi << 84 } while(greject < G4UniformRand()*grejsup); 83 // obtained for tsam=0 << 85 84 G4double gtmax = 2.0*(a1 + 1.0/ac); << 86 G4double sinteta = std::sqrt(sint2); 85 << 87 G4double phi = CLHEP::twopi*G4UniformRand(); 86 G4double tsam = 0.0; << 88 87 G4double gtr = 0.0; << 89 fLocalDirection.set(sinteta*std::cos(phi), sinteta*std::sin(phi), costeta); 88 << 89 //2) sampling. Eq. (2.31) of Penelope Manu << 90 // tsam = 1-std::cos(theta) << 91 // gtr = rejection function according to E << 92 do{ << 93 G4double rand = G4UniformRand(); << 94 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt << 95 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam) << 96 // Loop checking, 03-Aug-2015, Vladimir << 97 } while(G4UniformRand()*gtmax > gtr); << 98 << 99 costheta = 1.0 - tsam; << 100 << 101 G4double sint = std::sqrt(tsam*(2.0 - tsam << 102 G4double phi = CLHEP::twopi*G4UniformRand << 103 << 104 fLocalDirection.set(sint*std::cos(phi), si << 105 fLocalDirection.rotateUz(dp->GetMomentumDi 90 fLocalDirection.rotateUz(dp->GetMomentumDirection()); 106 } 91 } 107 return fLocalDirection; 92 return fLocalDirection; 108 } 93 } 109 94 110 void G4SauterGavrilaAngularDistribution::Print 95 void G4SauterGavrilaAngularDistribution::PrintGeneratorInformation() const 111 { 96 { 112 G4cout << "\n" << G4endl; 97 G4cout << "\n" << G4endl; 113 G4cout << "Non-polarized photoelectric effec 98 G4cout << "Non-polarized photoelectric effect angular generator." << G4endl; 114 G4cout << "The Sauter-Gavrila distribution f 99 G4cout << "The Sauter-Gavrila distribution for the K-shell is used."<<G4endl; 115 G4cout << "Originally developed by M.Maire f 100 G4cout << "Originally developed by M.Maire for Geant3" 116 << G4endl; << 101 << G4endl; 117 } 102 } 118 << 119 103