00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifndef G4MuMinusCaptureCascade_h
00037 #define G4MuMinusCaptureCascade_h 1
00038
00039 #include <CLHEP/Units/PhysicalConstants.h>
00040
00041 #include "globals.hh"
00042 #include "Randomize.hh"
00043 #include "G4ParticleDefinition.hh"
00044 #include "G4ThreeVector.hh"
00045
00046 class G4GHEKinematicsVector;
00047
00048 class G4MuMinusCaptureCascade
00049 {
00050 public:
00051
00052 G4MuMinusCaptureCascade();
00053
00054 ~G4MuMinusCaptureCascade();
00055
00056 G4int DoCascade(const G4double Z, const G4double A,
00057 G4GHEKinematicsVector* Cascade);
00058
00059 void DoBoundMuonMinusDecay(G4double Z,
00060 G4int* nCascade, G4GHEKinematicsVector* Cascade);
00061
00062 G4double GetKShellEnergy(G4double Z);
00063
00064 G4ThreeVector& GetRandomVec();
00065
00066 private:
00067
00068 G4double GetLinApprox(G4int N, const G4double* X, const G4double* Y,
00069 G4double Xuser);
00070
00071 void AddNewParticle(G4ParticleDefinition* aParticle,
00072 G4ThreeVector& Momentum,
00073 G4double mass,
00074 G4int* nParticle,
00075 G4GHEKinematicsVector* Cascade);
00076
00077
00078 G4MuMinusCaptureCascade& operator=(const G4MuMinusCaptureCascade &right);
00079 G4MuMinusCaptureCascade(const G4MuMinusCaptureCascade& );
00080
00081 G4double Emass, MuMass;
00082 G4ParticleDefinition* theElectron;
00083 G4ParticleDefinition* theGamma;
00084 G4ThreeVector randomVect;
00085 };
00086
00087
00088
00089 inline G4double G4MuMinusCaptureCascade::GetLinApprox(G4int N,
00090 const G4double* X,
00091 const G4double* Y,
00092 G4double Xuser)
00093 {
00094 G4double Yuser;
00095 if(Xuser <= X[0]) Yuser = Y[0];
00096 else if(Xuser >= X[N-1]) Yuser = Y[N-1];
00097 else {
00098 G4int i;
00099 for (i=1; i<N; i++){
00100 if(Xuser <= X[i]) break;
00101 }
00102
00103 if(Xuser == X[i]) Yuser = Y[i];
00104 else Yuser = Y[i-1] + (Y[i] - Y[i-1])*(Xuser - X[i-1])/(X[i] - X[i-1]);
00105 }
00106 return Yuser;
00107 }
00108
00109
00110
00111 inline G4ThreeVector& G4MuMinusCaptureCascade::GetRandomVec()
00112 {
00113
00114
00115
00116 G4double cost = 2.0 * G4UniformRand() - 1.0;
00117 G4double sint = std::sqrt((1.0 - cost)*(1.0 + cost));
00118 G4double Phi = CLHEP::twopi * G4UniformRand();
00119 randomVect = G4ThreeVector(sint * std::cos(Phi), sint * std::sin(Phi), cost);
00120 return randomVect;
00121 }
00122
00123 #endif
00124
00125