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 // 27 // 28 #ifndef G4FermiMomentum_h 29 #define G4FermiMomentum_h 1 30 31 #include "globals.hh" 32 #include "G4ThreeVector.hh" 33 #include "Randomize.hh" 34 #include "G4Pow.hh" 35 36 class G4FermiMomentum 37 { 38 39 public: 40 G4FermiMomentum(); 41 ~G4FermiMomentum(); 42 43 inline void Init(G4int anA, G4int aZ) {theA = anA; theZ = aZ;} 44 45 inline G4double GetFermiMomentum(G4double density) 46 { 47 return constofpmax * cbrt(density * theA); 48 } 49 50 inline G4ThreeVector GetMomentum(G4double density, 51 G4double maxMomentum=-1.) 52 { 53 if (maxMomentum < 0 ) maxMomentum=GetFermiMomentum(density); 54 G4ThreeVector p; 55 56 do { 57 p=G4ThreeVector(2.*G4UniformRand()-1., 58 2.*G4UniformRand()-1., 59 2.*G4UniformRand()-1.); 60 } while ( p.mag() > 1. ); /* Loop checking, 30-Oct-2015, G.Folger */ 61 return p*maxMomentum; 62 } 63 64 private: 65 66 G4double cbrt(G4double x) { return G4Pow::GetInstance()->A13(x); } 67 68 private: 69 70 G4int theA; 71 G4int theZ; 72 // pmax= hbar * c * ( 3* pi**2 * rho )**(1/3) = 73 // hbar * c * ( 3* pi**2 )**(1/3) * rho**(1/3)= 74 // constofpmax * rho**(1/3) 75 G4double constofpmax; 76 77 }; 78 79 #endif 80 81