ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronBetaDecayChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NeutronBetaDecayChannel.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 //
26 //
27 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class header file
31 //
32 // History: first implementation, based on object model of
33 // 18 Sep 2001 H.Kurashige
34 //---
35 // Fix energy of proton and neutrino May 2011 H.Kurashige
36 // Fix direction of proton and neutrino Nov 2013 H.Kurashige
37 // ------------------------------------------------------------
38 
39 #include "G4ParticleDefinition.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4DecayProducts.hh"
43 #include "G4VDecayChannel.hh"
45 #include "Randomize.hh"
46 #include "G4RotationMatrix.hh"
47 #include "G4LorentzVector.hh"
48 #include "G4LorentzRotation.hh"
49 
51  :G4VDecayChannel(),
52  aENuCorr(-0.102)
53 {
54 }
55 
57  const G4String& theParentName,
58  G4double theBR)
59  :G4VDecayChannel("Neutron Decay"),
60  aENuCorr(-0.102)
61 {
62  // set names for daughter particles
63  if (theParentName == "neutron") {
64  SetBR(theBR);
65  SetParent("neutron");
67  SetDaughter(0, "e-");
68  SetDaughter(1, "anti_nu_e");
69  SetDaughter(2, "proton");
70  } else if (theParentName == "anti_neutron") {
71  SetBR(theBR);
72  SetParent("anti_neutron");
74  SetDaughter(0, "e+");
75  SetDaughter(1, "nu_e");
76  SetDaughter(2, "anti_proton");
77  } else {
78 #ifdef G4VERBOSE
79  if (GetVerboseLevel()>0) {
80  G4cout << "G4NeutronBetaDecayChannel:: constructor :";
81  G4cout << " parent particle is not neutron but ";
82  G4cout << theParentName << G4endl;
83  }
84 #endif
85  }
86 }
87 
89 {
90 }
91 
93  : G4VDecayChannel(right),
94  aENuCorr(-0.102)
95 {
96 }
97 
98 
100 {
101  if (this != &right) {
103  verboseLevel = right.verboseLevel;
104  rbranch = right.rbranch;
105 
106  // copy parent name
107  parent_name = new G4String(*right.parent_name);
108 
109  // clear daughters_name array
111 
112  // recreate array
114  if ( numberOfDaughters >0 ) {
117  //copy daughters name
118  for (G4int index=0; index < numberOfDaughters; index++) {
119  daughters_name[index] = new G4String(*right.daughters_name[index]);
120  }
121  }
122  }
123  return *this;
124 }
125 
127 {
128  // This class describes free neutron beta decay kinemtics.
129  // This version neglects neutron/electron polarization
130  // without Coulomb effect
131 
132 #ifdef G4VERBOSE
133  if (GetVerboseLevel()>1) G4cout << "G4NeutronBetaDecayChannel::DecayIt ";
134 #endif
135 
138 
139  // parent mass
140  G4double parentmass = G4MT_parent->GetPDGMass();
141 
142  //daughters'mass
143  G4double daughtermass[3];
144  G4double sumofdaughtermass = 0.0;
145  for (G4int index=0; index<3; index++){
146  daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
147  sumofdaughtermass += daughtermass[index];
148  }
149  G4double xmax = parentmass-sumofdaughtermass;
150 
151  //create parent G4DynamicParticle at rest
152  G4ThreeVector dummy;
153  G4DynamicParticle * parentparticle = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
154 
155  //create G4Decayproducts
156  G4DecayProducts *products = new G4DecayProducts(*parentparticle);
157  delete parentparticle;
158 
159  // calculate daughter momentum
160  G4double daughtermomentum[3];
161 
162  // calcurate electron energy
163  G4double x; // Ee
164  G4double p; // Pe
165  G4double dm = daughtermass[0]; //Me
166  G4double w; // cosine of e-nu angle
167  G4double r;
168  G4double r0;
169  const size_t MAX_LOOP=10000;
170  for (size_t loop_counter=0; loop_counter <MAX_LOOP; ++loop_counter){
171  x = xmax*G4UniformRand();
172  p = std::sqrt(x*(x+2.0*dm));
173  w = 1.0-2.0*G4UniformRand();
174  r = p*(x+dm)*(xmax-x)*(xmax-x)*(1.0+aENuCorr*p/(x+dm)*w);
175  r0 = G4UniformRand()*(xmax+dm)*(xmax+dm)*xmax*xmax*(1.0+aENuCorr);
176  if (r > r0) break;
177  }
178 
179  //create daughter G4DynamicParticle
180  // rotation materix to lab frame
181  G4double costheta = 2.*G4UniformRand()-1.0;
182  G4double theta = std::acos(costheta)*rad;
184  G4RotationMatrix rm;
185  rm.rotateY(theta);
186  rm.rotateZ(phi);
187 
188  // daughter 0 (electron) in Z direction
189  daughtermomentum[0] = p;
190  G4ThreeVector direction0(0.0, 0.0, 1.0);
191  direction0 = rm * direction0;
192  G4DynamicParticle * daughterparticle0
193  = new G4DynamicParticle( G4MT_daughters[0], direction0*daughtermomentum[0]);
194  products->PushProducts(daughterparticle0);
195 
196  // daughter 1 (nutrino) in XZ plane
197  G4double eNu; // Enu
198  eNu = (parentmass-daughtermass[2])*(parentmass+daughtermass[2])+(dm*dm)-2.*parentmass*(x+dm);
199  eNu /= 2.*(parentmass+p*w-(x+dm));
200  G4double cosn = w;
201  G4double phin = twopi*G4UniformRand()*rad;
202  G4double sinn = std::sqrt((1.0-cosn)*(1.0+cosn));
203 
204  G4ThreeVector direction1(sinn*std::cos(phin), sinn*std::sin(phin), cosn);
205  direction1 = rm * direction1;
206  G4DynamicParticle * daughterparticle1
207  = new G4DynamicParticle( G4MT_daughters[1], direction1*eNu);
208  products->PushProducts(daughterparticle1);
209 
210  // daughter 2 (proton) at REST
211  G4double eP; // Eproton
212  eP = parentmass-eNu-(x+dm)-daughtermass[2];
213  G4double pPx = -eNu*sinn;
214  G4double pPz = -p-eNu*cosn;
215  G4double pP = std::sqrt(eP*(eP+2.*daughtermass[2]));
216  G4ThreeVector direction2(pPx/pP*std::cos(phin), pPx/pP*std::sin(phin), pPz/pP);
217  direction2 = rm * direction2;
218  G4DynamicParticle * daughterparticle2
219  = new G4DynamicParticle( G4MT_daughters[2], direction2*pP);
220  products->PushProducts(daughterparticle2);
221 
222 
223  // output message
224 #ifdef G4VERBOSE
225  if (GetVerboseLevel()>1) {
226  G4cout << "G4NeutronBetaDecayChannel::DecayIt ";
227  G4cout << " create decay products in rest frame " <<G4endl;
228  products->DumpInfo();
229  }
230 #endif
231  return products;
232 }
233 
234 
235 
236 
237 
238