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 // G4MuonRadiativeDecayChannelWithSpin 27 // 28 // Class description: 29 // 30 // This class describes radiative muon decay kinematics, but 31 // gives incorrect energy spectrum for neutrinos. 32 // Samples Radiative Muon Decay. 33 // References: 34 // - TRIUMF/TWIST Technote TN-55: 35 // "Radiative muon decay" by P. Depommier and A. Vacheret 36 // - Yoshitaka Kuno and Yasuhiro Okada 37 // "Muon Decays and Physics Beyond the Standard Model" 38 // Rev. Mod. Phys. 73, 151 (2001) 39 40 // Author: P.Gumplinger - Triumf, 25 July 2007 41 // Revision: D. Mingming - Center for HEP, Tsinghua Univ., 10 August 2011 42 // -------------------------------------------------------------------- 43 #ifndef G4MuonRadiativeDecayChannelWithSpin_hh 44 #define G4MuonRadiativeDecayChannelWithSpin_hh 1 45 46 #include "G4ThreeVector.hh" 47 #include "G4VDecayChannel.hh" 48 #include "Randomize.hh" 49 #include "globals.hh" 50 51 #include <CLHEP/Units/PhysicalConstants.h> 52 53 class G4MuonRadiativeDecayChannelWithSpin : public G4VDecayChannel 54 { 55 public: 56 G4MuonRadiativeDecayChannelWithSpin(const G4String& theParentName, G4double theBR); 57 ~G4MuonRadiativeDecayChannelWithSpin() override = default; 58 59 G4DecayProducts* DecayIt(G4double) override; 60 61 protected: 62 // Copy constructor and assignment operator 63 G4MuonRadiativeDecayChannelWithSpin(const G4MuonRadiativeDecayChannelWithSpin&) = default; 64 G4MuonRadiativeDecayChannelWithSpin& operator=(const G4MuonRadiativeDecayChannelWithSpin&); 65 66 private: 67 G4MuonRadiativeDecayChannelWithSpin() = default; 68 69 G4double fron(G4double Pmu, G4double x, G4double y, G4double cthetaE, G4double cthetaG, 70 G4double cthetaEG); 71 72 // Generates random vectors, uniformly distributed over the surface 73 // of a sphere of given radius 74 void rn3dim(G4double& x, G4double& y, G4double& z, G4double xlong); 75 76 G4double atan4(G4double x, G4double y); 77 }; 78 79 // ------------------------ 80 // Inline methods 81 // ------------------------ 82 83 inline void G4MuonRadiativeDecayChannelWithSpin::rn3dim(G4double& x, G4double& y, G4double& z, 84 G4double xlong) 85 { 86 G4double a = 0.; 87 G4double b = 0.; 88 G4double c = 0.; 89 G4double r = 0.; 90 91 do { 92 a = G4UniformRand() - 0.5; 93 b = G4UniformRand() - 0.5; 94 c = G4UniformRand() - 0.5; 95 r = a * a + b * b + c * c; 96 } while (r > 0.25); // Loop checking, 09.08.2015, K.Kurashige 97 98 G4double rinv = xlong / (std::sqrt(r)); 99 x = a * rinv; 100 y = b * rinv; 101 z = c * rinv; 102 103 return; 104 } 105 106 inline G4double G4MuonRadiativeDecayChannelWithSpin::atan4(G4double x, G4double y) 107 { 108 G4double phi = 0.; 109 110 if (x == 0. && y > 0.) { 111 phi = 0.5 * CLHEP::pi; 112 } 113 else if (x == 0. && y < 0.) { 114 phi = 1.5 * CLHEP::pi; 115 } 116 else if (y == 0. && x > 0.) { 117 phi = 0.; 118 } 119 else if (y == 0. && x < 0.) { 120 phi = CLHEP::pi; 121 } 122 else if (x > 0.) { 123 phi = std::atan(y / x); 124 } 125 else if (x < 0.) { 126 phi = std::atan(y / x) + CLHEP::pi; 127 } 128 129 return phi; 130 } 131 132 #endif 133