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