ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NeutronDecay.cc
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 //
27 // //
28 // File: G4NeutronDecay.cc //
29 // Author: L.G. Sarmiento (Lund) //
30 // Date: 10 October 2015 //
31 // //
33 
34 #include "G4NeutronDecay.hh"
35 #include "G4IonTable.hh"
36 #include "Randomize.hh"
37 #include "G4ThreeVector.hh"
38 #include "G4DynamicParticle.hh"
39 #include "G4DecayProducts.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include <iostream>
43 #include <iomanip>
44 
46  const G4double& branch, const G4double& Qvalue,
47  const G4double& excitationE,
48  const G4Ions::G4FloatLevelBase& flb)
49  : G4NuclearDecay("neutron decay", Neutron, excitationE, flb), transitionQ(Qvalue)
50 {
51  SetParent(theParentNucleus); // Store name of parent nucleus, delete G4MT_parent
52  SetBR(branch);
53 
55  G4IonTable* theIonTable =
57  G4int daughterZ = theParentNucleus->GetAtomicNumber();
58  G4int daughterA = theParentNucleus->GetAtomicMass() - 1;
59  SetDaughter(0, theIonTable->GetIon(daughterZ, daughterA, excitationE, flb) );
60  SetDaughter(1, "neutron");
61 }
62 
63 
65 {}
66 
67 
69 {
70  // Fill G4MT_parent with theParentNucleus (stored by SetParent in ctor)
72 
73  // Fill G4MT_daughters with neutron and residual nucleus (stored by SetDaughter)
75 
76  G4double neutronMass = G4MT_daughters[1]->GetPDGMass();
77  // Excitation energy included in PDG mass
78  G4double nucleusMass = G4MT_daughters[0]->GetPDGMass();
79 
80  // Q value was calculated from atomic masses.
81  // Use it to get correct neutron energy.
82  G4double cmMomentum = std::sqrt(transitionQ*(transitionQ + 2.*neutronMass)*
83  (transitionQ + 2.*nucleusMass)*
84  (transitionQ + 2.*neutronMass + 2.*nucleusMass) )/
85  (transitionQ + neutronMass + nucleusMass)/2.;
86 
87  // Set up final state
88  // parentParticle is set at rest here because boost with correct momentum
89  // is done later
90  G4DynamicParticle parentParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
91  G4DecayProducts* products = new G4DecayProducts(parentParticle);
92 
93  G4double costheta = 2.*G4UniformRand()-1.0;
94  G4double sintheta = std::sqrt(1.0 - costheta*costheta);
96  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),
97  costheta);
98 
99  G4double KE = std::sqrt(cmMomentum*cmMomentum + neutronMass*neutronMass)
100  - neutronMass;
101  G4DynamicParticle* daughterparticle =
102  new G4DynamicParticle(G4MT_daughters[1], direction, KE, neutronMass);
103  products->PushProducts(daughterparticle);
104 
105  KE = std::sqrt(cmMomentum*cmMomentum + nucleusMass*nucleusMass) - nucleusMass;
106  daughterparticle =
107  new G4DynamicParticle(G4MT_daughters[0], -1.0*direction, KE, nucleusMass);
108  products->PushProducts(daughterparticle);
109 
110  // Energy conservation check
111  // For neutron decays, do final energy check against reaction Q value
112  // which is well-measured using atomic mass differences. Nuclear masses
113  // should not be used since they are not usually directly measured and we
114  // always decay atoms and not fully stripped nuclei.
115  /*
116  G4int nProd = products->entries();
117  G4DynamicParticle* temp = 0;
118  G4double Esum = 0.0;
119  for (G4int i = 0; i < nProd; i++) {
120  temp = products->operator[](i);
121  Esum += temp->GetKineticEnergy();
122  }
123  G4double eCons = (transitionQ - Esum)/keV;
124  if (eCons > 1.e-07) G4cout << " Neutron decay check: Ediff (keV) = " << eCons << G4endl;
125  */
126  return products;
127 }
128 
129 
131 {
132  G4cout << " G4NeutronDecay for parent nucleus " << GetParentName() << G4endl;
133  G4cout << " decays to " << GetDaughterName(0) << " + " << GetDaughterName(1)
134  << " with branching ratio " << GetBR() << "% and Q value "
135  << transitionQ << G4endl;
136 }
137