| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // G4PionRadiativeDecayChannel class implement
27 // GEANT 4 class header file
28 //
29 // Author: P.Gumplinger, 30 July 2007
30 // Reference: M. Blecher, TRIUMF/PIENU Technot
31 // "Inclusion of pi->enug in the Mo
32 // -------------------------------------------
33
34 #include "G4PionRadiativeDecayChannel.hh"
35
36 #include "G4DecayProducts.hh"
37 #include "G4LorentzVector.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "Randomize.hh"
41
42 namespace
43 {
44 const G4double beta = 3.6612e-03;
45 const G4double cib = 1.16141e-03;
46 const G4double csdp = 3.45055e-02;
47 const G4double csdm = 5.14122e-03;
48 const G4double cif = 4.63543e-05;
49 const G4double cig = 1.78928e-05;
50 const G4double xl = 2. * 0.1 * MeV / 139.57 *
51 const G4double yl = ((1. - xl) + std::sqrt((1
52
53 const G4double xu = 1. - (yl - std::sqrt(yl *
54 const G4double yu = 1. + beta * beta;
55
56 inline G4double D2W(const G4double x, const G4
57 {
58 return cib * (1. - y) * (1. + ((1. - x) * (1
59 + csdp * (1. - x) * ((x + y - 1.) * (
60 + cif * (x - 1.) * (1. - y) / x + cig
61 }
62
63 const G4double d2wmax = D2W(xl, yl);
64 } // namespace
65
66 G4PionRadiativeDecayChannel::G4PionRadiativeDe
67
68 : G4VDecayChannel("Radiative Pion Decay", 1)
69 {
70 // set names for daughter particles
71 if (theParentName == "pi+") {
72 SetBR(theBR);
73 SetParent("pi+");
74 SetNumberOfDaughters(3);
75 SetDaughter(0, "e+");
76 SetDaughter(1, "gamma");
77 SetDaughter(2, "nu_e");
78 }
79 else if (theParentName == "pi-") {
80 SetBR(theBR);
81 SetParent("pi-");
82 SetNumberOfDaughters(3);
83 SetDaughter(0, "e-");
84 SetDaughter(1, "gamma");
85 SetDaughter(2, "anti_nu_e");
86 }
87 else {
88 #ifdef G4VERBOSE
89 if (GetVerboseLevel() > 0) {
90 G4cout << "G4RadiativePionDecayChannel::
91 G4cout << "Parent particle is not charge
92 G4cout << theParentName << G4endl;
93 }
94 #endif
95 }
96 }
97
98 G4PionRadiativeDecayChannel&
99 G4PionRadiativeDecayChannel::operator=(const G
100 {
101 if (this != &right) {
102 kinematics_name = right.kinematics_name;
103 verboseLevel = right.verboseLevel;
104 rbranch = right.rbranch;
105
106 // copy parent name
107 parent_name = new G4String(*right.parent_n
108
109 // clear daughters_name array
110 ClearDaughtersName();
111
112 // recreate array
113 numberOfDaughters = right.numberOfDaughter
114 if (numberOfDaughters > 0) {
115 if (daughters_name != nullptr) ClearDaug
116 daughters_name = new G4String*[numberOfD
117 // copy daughters name
118 for (G4int index = 0; index < numberOfDa
119 daughters_name[index] = new G4String(*
120 }
121 }
122 }
123 return *this;
124 }
125
126 G4DecayProducts* G4PionRadiativeDecayChannel::
127 {
128 #ifdef G4VERBOSE
129 if (GetVerboseLevel() > 1) G4cout << "G4Pion
130 #endif
131
132 CheckAndFillParent();
133 CheckAndFillDaughters();
134
135 // parent mass
136 G4double parentmass = G4MT_parent->GetPDGMas
137
138 G4double EMPI = parentmass;
139
140 // daughters'mass
141 const G4int N_DAUGHTER = 3;
142 G4double daughtermass[N_DAUGHTER];
143 // G4double sumofdaughtermass = 0.0;
144 for (G4int index = 0; index < N_DAUGHTER; ++
145 daughtermass[index] = G4MT_daughters[index
146 // sumofdaughtermass += daughtermass[index
147 }
148
149 G4double EMASS = daughtermass[0];
150
151 // create parent G4DynamicParticle at rest
152 G4ThreeVector dummy;
153 auto parentparticle = new G4DynamicParticle(
154 // create G4Decayproducts
155 auto products = new G4DecayProducts(*parentp
156 delete parentparticle;
157
158 G4double x, y;
159
160 const std::size_t MAX_LOOP = 1000;
161
162 for (std::size_t loop_counter1 = 0; loop_cou
163 for (std::size_t loop_counter2 = 0; loop_c
164 x = xl + G4UniformRand() * (xu - xl);
165 y = yl + G4UniformRand() * (yu - yl);
166 if (x + y > 1.) break;
167 }
168 G4double d2w = D2W(x, y);
169 if (d2w > G4UniformRand() * d2wmax) break;
170 }
171
172 // Calculate the angle between positron and
173 //
174 G4double cthetaGE =
175 (y * (x - 2.) + 2. * (1. - x + beta * beta
176
177 G4double G = x * EMPI / 2.;
178 G4double E = y * EMPI / 2.;
179
180 if (E < EMASS) E = EMASS;
181
182 // calculate daughter momentum
183 G4double daughtermomentum[2];
184
185 daughtermomentum[0] = std::sqrt(E * E - EMAS
186
187 G4double cthetaE = 2. * G4UniformRand() - 1.
188 G4double sthetaE = std::sqrt(1. - cthetaE *
189
190 G4double phiE = twopi * G4UniformRand() * ra
191 G4double cphiE = std::cos(phiE);
192 G4double sphiE = std::sin(phiE);
193
194 // Coordinates of the decay positron
195 //
196 G4double px = sthetaE * cphiE;
197 G4double py = sthetaE * sphiE;
198 G4double pz = cthetaE;
199
200 G4ThreeVector direction0(px, py, pz);
201
202 auto daughterparticle0 =
203 new G4DynamicParticle(G4MT_daughters[0], d
204
205 products->PushProducts(daughterparticle0);
206
207 daughtermomentum[1] = G;
208
209 G4double sthetaGE = std::sqrt(1. - cthetaGE
210
211 G4double phiGE = twopi * G4UniformRand() * r
212 G4double cphiGE = std::cos(phiGE);
213 G4double sphiGE = std::sin(phiGE);
214
215 // Coordinates of the decay gamma with respe
216 //
217 px = sthetaGE * cphiGE;
218 py = sthetaGE * sphiGE;
219 pz = cthetaGE;
220
221 G4ThreeVector direction1(px, py, pz);
222
223 direction1.rotateUz(direction0);
224
225 auto daughterparticle1 =
226 new G4DynamicParticle(G4MT_daughters[1], d
227
228 products->PushProducts(daughterparticle1);
229
230 // output message
231 #ifdef G4VERBOSE
232 if (GetVerboseLevel() > 1) {
233 G4cout << "G4PionRadiativeDecayChannel::De
234 G4cout << " create decay products in rest
235 products->DumpInfo();
236 }
237 #endif
238
239 return products;
240 }
241