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 // G4NeutronBetaDecayChannel class implementation 27 // 28 // Author: H.Kurashige, 18 September 2001 29 // -------------------------------------------------------------------- 30 31 #include "G4NeutronBetaDecayChannel.hh" 32 33 #include "G4DecayProducts.hh" 34 #include "G4LorentzRotation.hh" 35 #include "G4LorentzVector.hh" 36 #include "G4ParticleDefinition.hh" 37 #include "G4PhysicalConstants.hh" 38 #include "G4RotationMatrix.hh" 39 #include "G4SystemOfUnits.hh" 40 #include "G4VDecayChannel.hh" 41 #include "Randomize.hh" 42 43 G4NeutronBetaDecayChannel::G4NeutronBetaDecayChannel(const G4String& theParentName, G4double theBR) 44 : G4VDecayChannel("Neutron Decay") 45 { 46 // set names for daughter particles 47 if (theParentName == "neutron") { 48 SetBR(theBR); 49 SetParent("neutron"); 50 SetNumberOfDaughters(3); 51 SetDaughter(0, "e-"); 52 SetDaughter(1, "anti_nu_e"); 53 SetDaughter(2, "proton"); 54 } 55 else if (theParentName == "anti_neutron") { 56 SetBR(theBR); 57 SetParent("anti_neutron"); 58 SetNumberOfDaughters(3); 59 SetDaughter(0, "e+"); 60 SetDaughter(1, "nu_e"); 61 SetDaughter(2, "anti_proton"); 62 } 63 else { 64 #ifdef G4VERBOSE 65 if (GetVerboseLevel() > 0) { 66 G4cout << "G4NeutronBetaDecayChannel:: constructor :"; 67 G4cout << " parent particle is not neutron but "; 68 G4cout << theParentName << G4endl; 69 } 70 #endif 71 } 72 } 73 74 G4NeutronBetaDecayChannel::G4NeutronBetaDecayChannel(const G4NeutronBetaDecayChannel& right) 75 : G4VDecayChannel(right) 76 {} 77 78 G4NeutronBetaDecayChannel& 79 G4NeutronBetaDecayChannel::operator=(const G4NeutronBetaDecayChannel& right) 80 { 81 if (this != &right) { 82 kinematics_name = right.kinematics_name; 83 verboseLevel = right.verboseLevel; 84 rbranch = right.rbranch; 85 86 // copy parent name 87 delete parent_name; 88 parent_name = new G4String(*right.parent_name); 89 90 // clear daughters_name array 91 ClearDaughtersName(); 92 93 // recreate array 94 numberOfDaughters = right.numberOfDaughters; 95 if (numberOfDaughters > 0) { 96 daughters_name = new G4String*[numberOfDaughters]; 97 // copy daughters name 98 for (G4int index = 0; index < numberOfDaughters; ++index) { 99 daughters_name[index] = new G4String(*right.daughters_name[index]); 100 } 101 } 102 } 103 return *this; 104 } 105 106 G4DecayProducts* G4NeutronBetaDecayChannel::DecayIt(G4double) 107 { 108 // This class describes free neutron beta decay kinematics. 109 // This version neglects neutron/electron polarization 110 // without Coulomb effect 111 112 #ifdef G4VERBOSE 113 if (GetVerboseLevel() > 1) G4cout << "G4NeutronBetaDecayChannel::DecayIt "; 114 #endif 115 116 CheckAndFillParent(); 117 CheckAndFillDaughters(); 118 119 // parent mass 120 G4double parentmass = G4MT_parent->GetPDGMass(); 121 122 // daughters'mass 123 G4double daughtermass[3]; 124 G4double sumofdaughtermass = 0.0; 125 for (G4int index = 0; index < 3; ++index) { 126 daughtermass[index] = G4MT_daughters[index]->GetPDGMass(); 127 sumofdaughtermass += daughtermass[index]; 128 } 129 G4double xmax = parentmass - sumofdaughtermass; 130 131 // create parent G4DynamicParticle at rest 132 G4ThreeVector dummy; 133 auto parentparticle = new G4DynamicParticle(G4MT_parent, dummy, 0.0); 134 135 // create G4Decayproducts 136 auto products = new G4DecayProducts(*parentparticle); 137 delete parentparticle; 138 139 // calculate daughter momentum 140 G4double daughtermomentum[3]; 141 142 // calcurate electron energy 143 G4double x; // Ee 144 G4double p; // Pe 145 G4double dm = daughtermass[0]; // Me 146 G4double w; // cosine of e-nu angle 147 G4double r; 148 G4double r0; 149 const std::size_t MAX_LOOP = 10000; 150 for (std::size_t loop_counter = 0; loop_counter < MAX_LOOP; ++loop_counter) { 151 x = xmax * G4UniformRand(); 152 p = std::sqrt(x * (x + 2.0 * dm)); 153 w = 1.0 - 2.0 * G4UniformRand(); 154 r = p * (x + dm) * (xmax - x) * (xmax - x) * (1.0 + aENuCorr * p / (x + dm) * w); 155 r0 = G4UniformRand() * (xmax + dm) * (xmax + dm) * xmax * xmax * (1.0 + aENuCorr); 156 if (r > r0) break; 157 } 158 159 // create daughter G4DynamicParticle 160 // rotation materix to lab frame 161 // 162 G4double costheta = 2. * G4UniformRand() - 1.0; 163 G4double theta = std::acos(costheta) * rad; 164 G4double phi = twopi * G4UniformRand() * rad; 165 G4RotationMatrix rm; 166 rm.rotateY(theta); 167 rm.rotateZ(phi); 168 169 // daughter 0 (electron) in Z direction 170 daughtermomentum[0] = p; 171 G4ThreeVector direction0(0.0, 0.0, 1.0); 172 direction0 = rm * direction0; 173 auto daughterparticle0 = 174 new G4DynamicParticle(G4MT_daughters[0], direction0 * daughtermomentum[0]); 175 products->PushProducts(daughterparticle0); 176 177 // daughter 1 (nutrino) in XZ plane 178 G4double eNu; // Enu 179 eNu = (parentmass - daughtermass[2]) * (parentmass + daughtermass[2]) + (dm * dm) 180 - 2. * parentmass * (x + dm); 181 eNu /= 2. * (parentmass + p * w - (x + dm)); 182 G4double cosn = w; 183 G4double phin = twopi * G4UniformRand() * rad; 184 G4double sinn = std::sqrt((1.0 - cosn) * (1.0 + cosn)); 185 186 G4ThreeVector direction1(sinn * std::cos(phin), sinn * std::sin(phin), cosn); 187 direction1 = rm * direction1; 188 auto daughterparticle1 = new G4DynamicParticle(G4MT_daughters[1], direction1 * eNu); 189 products->PushProducts(daughterparticle1); 190 191 // daughter 2 (proton) at REST 192 G4double eP; // Eproton 193 eP = parentmass - eNu - (x + dm) - daughtermass[2]; 194 G4double pPx = -eNu * sinn; 195 G4double pPz = -p - eNu * cosn; 196 G4double pP = std::sqrt(eP * (eP + 2. * daughtermass[2])); 197 G4ThreeVector direction2(pPx / pP * std::cos(phin), pPx / pP * std::sin(phin), pPz / pP); 198 direction2 = rm * direction2; 199 auto daughterparticle2 = new G4DynamicParticle(G4MT_daughters[2], direction2 * pP); 200 products->PushProducts(daughterparticle2); 201 202 // output message 203 #ifdef G4VERBOSE 204 if (GetVerboseLevel() > 1) { 205 G4cout << "G4NeutronBetaDecayChannel::DecayIt "; 206 G4cout << " create decay products in rest frame " << G4endl; 207 products->DumpInfo(); 208 } 209 #endif 210 return products; 211 } 212