| 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 // G4MuonRadiativeDecayChannelWithSpin class i << 26 // ------------------------------------------------------------
>> 27 // GEANT 4 class header file
>> 28 //
>> 29 // History:
>> 30 // 01 August 2007 P.Gumplinger
>> 31 // 10 August 2011 D. Mingming - Center for HEP, Tsinghua Univ.
>> 32 // References:
>> 33 // TRIUMF/TWIST Technote TN-55:
>> 34 // "Radiative muon decay" by P. Depommier and A. Vacheret
>> 35 // ------------------------------------------------------
>> 36 // Yoshitaka Kuno and Yasuhiro Okada
>> 37 // "Muon Decays and Physics Beyond the Standard Model"
>> 38 // Rev. Mod. Phys. 73, 151 (2001)
>> 39 //
>> 40 // ------------------------------------------------------------
>> 41 //
>> 42 //
27 // 43 //
28 // References: <<
29 // - TRIUMF/TWIST Technote TN-55: <<
30 // "Radiative muon decay" by P. Depommier an <<
31 // - Yoshitaka Kuno and Yasuhiro Okada <<
32 // "Muon Decays and Physics Beyond the Stand <<
33 // Rev. Mod. Phys. 73, 151 (2001) <<
34 // <<
35 // Author: P.Gumplinger - Triumf, 25 July 2007 <<
36 // Revision: D.Mingming - Center for HEP, Tsin <<
37 // ------------------------------------------- <<
38 44
39 #include "G4MuonRadiativeDecayChannelWithSpin. 45 #include "G4MuonRadiativeDecayChannelWithSpin.hh"
40 46
41 #include "G4DecayProducts.hh" <<
42 #include "G4LorentzVector.hh" <<
43 #include "G4PhysicalConstants.hh" 47 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh" 48 #include "G4SystemOfUnits.hh"
45 #include "Randomize.hh" 49 #include "Randomize.hh"
>> 50 #include "G4DecayProducts.hh"
>> 51 #include "G4LorentzVector.hh"
46 52
47 G4MuonRadiativeDecayChannelWithSpin::G4MuonRad << 53 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin()
48 const G4String& theParentName, G4double theB << 54 : G4VDecayChannel()
49 : G4VDecayChannel("Radiative Muon Decay", 1) << 55 {
>> 56 }
>> 57
>> 58 G4MuonRadiativeDecayChannelWithSpin::
>> 59 G4MuonRadiativeDecayChannelWithSpin(const G4String& theParentName,
>> 60 G4double theBR)
>> 61 : G4VDecayChannel("Radiative Muon Decay",1)
50 { 62 {
51 // set names for daughter particles 63 // set names for daughter particles
52 if (theParentName == "mu+") { 64 if (theParentName == "mu+") {
53 SetBR(theBR); 65 SetBR(theBR);
54 SetParent("mu+"); 66 SetParent("mu+");
55 SetNumberOfDaughters(4); 67 SetNumberOfDaughters(4);
56 SetDaughter(0, "e+"); 68 SetDaughter(0, "e+");
57 SetDaughter(1, "gamma"); 69 SetDaughter(1, "gamma");
58 SetDaughter(2, "nu_e"); 70 SetDaughter(2, "nu_e");
59 SetDaughter(3, "anti_nu_mu"); 71 SetDaughter(3, "anti_nu_mu");
60 } << 72 } else if (theParentName == "mu-") {
61 else if (theParentName == "mu-") { <<
62 SetBR(theBR); 73 SetBR(theBR);
63 SetParent("mu-"); 74 SetParent("mu-");
64 SetNumberOfDaughters(4); 75 SetNumberOfDaughters(4);
65 SetDaughter(0, "e-"); 76 SetDaughter(0, "e-");
66 SetDaughter(1, "gamma"); 77 SetDaughter(1, "gamma");
67 SetDaughter(2, "anti_nu_e"); 78 SetDaughter(2, "anti_nu_e");
68 SetDaughter(3, "nu_mu"); 79 SetDaughter(3, "nu_mu");
69 } << 80 } else {
70 else { <<
71 #ifdef G4VERBOSE 81 #ifdef G4VERBOSE
72 if (GetVerboseLevel() > 0) { << 82 if (GetVerboseLevel()>0) {
73 G4cout << "G4RadiativeMuonDecayChannel:: << 83 G4cout << "G4RadiativeMuonDecayChannel:: constructor :";
74 G4cout << " parent particle is not muon 84 G4cout << " parent particle is not muon but ";
75 G4cout << theParentName << G4endl; 85 G4cout << theParentName << G4endl;
76 } 86 }
77 #endif 87 #endif
78 } 88 }
79 } 89 }
80 90
81 G4MuonRadiativeDecayChannelWithSpin& << 91 G4MuonRadiativeDecayChannelWithSpin::~G4MuonRadiativeDecayChannelWithSpin()
82 G4MuonRadiativeDecayChannelWithSpin::operator= << 92 {
>> 93 }
>> 94
>> 95 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin(const G4MuonRadiativeDecayChannelWithSpin &right):
>> 96 G4VDecayChannel(right)
>> 97 {
>> 98 }
>> 99
>> 100 G4MuonRadiativeDecayChannelWithSpin & G4MuonRadiativeDecayChannelWithSpin::operator=(const G4MuonRadiativeDecayChannelWithSpin & right)
83 { 101 {
84 if (this != &right) { << 102 if (this != &right) {
85 kinematics_name = right.kinematics_name; 103 kinematics_name = right.kinematics_name;
86 verboseLevel = right.verboseLevel; 104 verboseLevel = right.verboseLevel;
87 rbranch = right.rbranch; 105 rbranch = right.rbranch;
88 106
89 // copy parent name 107 // copy parent name
90 parent_name = new G4String(*right.parent_n 108 parent_name = new G4String(*right.parent_name);
91 109
92 // clear daughters_name array 110 // clear daughters_name array
93 ClearDaughtersName(); 111 ClearDaughtersName();
94 112
95 // recreate array 113 // recreate array
96 numberOfDaughters = right.numberOfDaughter 114 numberOfDaughters = right.numberOfDaughters;
97 if (numberOfDaughters > 0) { << 115 if ( numberOfDaughters >0 ) {
98 if (daughters_name != nullptr) ClearDaug << 116 if (daughters_name !=0) ClearDaughtersName();
99 daughters_name = new G4String*[numberOfD 117 daughters_name = new G4String*[numberOfDaughters];
100 // copy daughters name << 118 //copy daughters name
101 for (G4int index = 0; index < numberOfDa << 119 for (G4int index=0; index < numberOfDaughters; index++) {
102 daughters_name[index] = new G4String(* << 120 daughters_name[index] = new G4String(*right.daughters_name[index]);
103 } 121 }
104 } 122 }
105 parent_polarization = right.parent_polariz 123 parent_polarization = right.parent_polarization;
106 } 124 }
107 return *this; 125 return *this;
108 } 126 }
109 127
110 G4DecayProducts* G4MuonRadiativeDecayChannelWi << 128
>> 129 G4DecayProducts *G4MuonRadiativeDecayChannelWithSpin::DecayIt(G4double)
111 { 130 {
>> 131
112 #ifdef G4VERBOSE 132 #ifdef G4VERBOSE
113 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 133 if (GetVerboseLevel()>1)
>> 134 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
114 #endif 135 #endif
115 136
116 CheckAndFillParent(); 137 CheckAndFillParent();
117 CheckAndFillDaughters(); 138 CheckAndFillDaughters();
118 139
119 // parent mass 140 // parent mass
120 G4double parentmass = G4MT_parent->GetPDGMas 141 G4double parentmass = G4MT_parent->GetPDGMass();
121 142
122 G4double EMMU = parentmass; 143 G4double EMMU = parentmass;
123 144
124 // daughters'mass << 145 //daughters'mass
125 G4double daughtermass[4]; << 146 G4double daughtermass[4];
126 // G4double sumofdaughtermass = 0.0; << 147 G4double sumofdaughtermass = 0.0;
127 for (G4int index = 0; index < 4; ++index) { << 148 for (G4int index=0; index<4; index++){
128 daughtermass[index] = G4MT_daughters[index 149 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
129 // sumofdaughtermass += daughtermass[index << 150 sumofdaughtermass += daughtermass[index];
130 } 151 }
131 152
132 G4double EMASS = daughtermass[0]; 153 G4double EMASS = daughtermass[0];
133 154
134 // create parent G4DynamicParticle at rest << 155 //create parent G4DynamicParticle at rest
135 G4ThreeVector dummy; 156 G4ThreeVector dummy;
136 auto parentparticle = new G4DynamicParticle( << 157 G4DynamicParticle * parentparticle =
137 << 158 new G4DynamicParticle( G4MT_parent, dummy, 0.0);
138 // create G4Decayproducts << 159 //create G4Decayproducts
139 auto products = new G4DecayProducts(*parentp << 160 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
140 delete parentparticle; 161 delete parentparticle;
141 162
142 G4double eps = EMASS / EMMU; << 163 G4int i = 0;
>> 164
>> 165 G4double eps = EMASS/EMMU;
143 166
144 G4double som0, x, y, xx, yy, zz; 167 G4double som0, x, y, xx, yy, zz;
145 G4double cthetaE, cthetaG, cthetaGE, phiE, p 168 G4double cthetaE, cthetaG, cthetaGE, phiE, phiG;
146 const std::size_t MAX_LOOP = 10000; << 169 const size_t MAX_LOOP=10000;
147 170
148 for (std::size_t loop_counter1 = 0; loop_cou << 171 for (size_t loop_counter1=0; loop_counter1 <MAX_LOOP; ++loop_counter1){
149 for (std::size_t loop_counter2 = 0; loop_c <<
150 // ------------------------------------- <<
151 // Build two vectors of random length an <<
152 // positron and the photon. <<
153 // x/y is the length of the vector, xx, <<
154 // phi is the azimutal angle, theta the <<
155 // ------------------------------------- <<
156 <<
157 // For the positron <<
158 // <<
159 x = G4UniformRand(); <<
160 172
161 rn3dim(xx, yy, zz, x); << 173 // leap1:
162 174
163 if (std::fabs((xx * xx) + (yy * yy) + (z << 175 i++;
164 G4cout << "Norm of x not correct" << G <<
165 } <<
166 176
167 phiE = atan4(xx, yy); << 177 // leap2:
168 cthetaE = zz / x; <<
169 G4double sthetaE = std::sqrt((xx * xx) + <<
170 <<
171 // What you get: <<
172 // <<
173 // x = positron energy <<
174 // phiE = azimutal angle of positron <<
175 // cthetaE = cosine of polar angle of po <<
176 // sthetaE = sine of polar angle of posi <<
177 // <<
178 //// G4cout << " x, xx, yy, zz " << x < <<
179 //// << yy < <<
180 //// G4cout << " phiE, cthetaE, sthetaE <<
181 //// <<
182 //// <<
183 <<
184 // For the photon <<
185 // <<
186 y = G4UniformRand(); <<
187 178
188 rn3dim(xx, yy, zz, y); << 179 for (size_t loop_counter2=0; loop_counter2 <MAX_LOOP; ++loop_counter2){
>> 180 //
>> 181 //--------------------------------------------------------------------------
>> 182 // Build two vectors of random length and random direction, for the
>> 183 // positron and the photon.
>> 184 // x/y is the length of the vector, xx, yy and zz the components,
>> 185 // phi is the azimutal angle, theta the polar angle.
>> 186 //--------------------------------------------------------------------------
>> 187 //
>> 188 // For the positron
>> 189 //
>> 190 x = G4UniformRand();
189 191
190 if (std::fabs((xx * xx) + (yy * yy) + (z << 192 rn3dim(xx,yy,zz,x);
191 G4cout << " Norm of y not correct " << <<
192 } <<
193 193
194 phiG = atan4(xx, yy); << 194 if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(x*x))>0.001){
195 cthetaG = zz / y; << 195 G4cout << "Norm of x not correct" << G4endl;
196 G4double sthetaG = std::sqrt((xx * xx) + << 196 }
197 <<
198 // What you get: <<
199 // <<
200 // y = photon energy <<
201 // phiG = azimutal angle of photon mo <<
202 // cthetaG = cosine of polar angle of ph <<
203 // sthetaG = sine of polar angle of phot <<
204 // <<
205 //// G4cout << " y, xx, yy, zz " << y < <<
206 //// << yy < <<
207 //// G4cout << " phiG, cthetaG, sthetaG <<
208 //// <<
209 //// <<
210 <<
211 // Calculate the angle between posi <<
212 // <<
213 cthetaGE = cthetaE * cthetaG + sthetaE * <<
214 <<
215 //// G4cout << x << " " << cthetaE << " <<
216 //// << y << " " << cthetaG << " <<
217 //// << cthetaGE <<
218 <<
219 G4double term0 = eps * eps; <<
220 G4double term1 = x * ((1.0 - eps) * (1.0 <<
221 G4double beta = <<
222 std::sqrt(x * ((1.0 - eps) * (1.0 - ep <<
223 / term1; <<
224 G4double delta = 1.0 - beta * cthetaGE; <<
225 <<
226 G4double term3 = y * (1.0 - (eps * eps)) <<
227 G4double term6 = term1 * delta * term3; <<
228 <<
229 G4double Qsqr = (1.0 - term1 - term3 + t <<
230 <<
231 // Check the kinematics. <<
232 // <<
233 if (Qsqr >= 0.0 && Qsqr <= 1.0) break; <<
234 <<
235 } // end loop count <<
236 <<
237 // Do the calculation for -1 muon polariza <<
238 // <<
239 G4double Pmu = -1.0; <<
240 if (GetParentName() == "mu-") { <<
241 Pmu = +1.0; <<
242 } <<
243 197
244 som0 = fron(Pmu, x, y, cthetaE, cthetaG, c << 198 phiE = atan4(xx,yy);
>> 199 cthetaE = zz/x;
>> 200 G4double sthetaE = std::sqrt((xx*xx)+(yy*yy))/x;
>> 201 //
>> 202 // What you get:
>> 203 //
>> 204 // x = positron energy
>> 205 // phiE = azimutal angle of positron momentum
>> 206 // cthetaE = cosine of polar angle of positron momentum
>> 207 // sthetaE = sine of polar angle of positron momentum
>> 208 //
>> 209 //// G4cout << " x, xx, yy, zz " << x << " " << xx << " "
>> 210 //// << yy << " " << zz << G4endl;
>> 211 //// G4cout << " phiE, cthetaE, sthetaE " << phiE << " "
>> 212 //// << cthetaE << " "
>> 213 //// << sthetaE << " " << G4endl;
>> 214 //
>> 215 //-----------------------------------------------------------------------
>> 216 //
>> 217 // For the photon
>> 218 //
>> 219 y = G4UniformRand();
>> 220
>> 221 rn3dim(xx,yy,zz,y);
245 222
246 // Sample the decay rate << 223 if(std::fabs((xx*xx)+(yy*yy)+(zz*zz)-(y*y))>0.001){
247 // << 224 G4cout << " Norm of y not correct " << G4endl;
248 if (G4UniformRand() * 177.0 <= som0) break << 225 }
>> 226
>> 227 phiG = atan4(xx,yy);
>> 228 cthetaG = zz/y;
>> 229 G4double sthetaG = std::sqrt((xx*xx)+(yy*yy))/y;
>> 230 //
>> 231 // What you get:
>> 232 //
>> 233 // y = photon energy
>> 234 // phiG = azimutal angle of photon momentum
>> 235 // cthetaG = cosine of polar angle of photon momentum
>> 236 // sthetaG = sine of polar angle of photon momentum
>> 237 //
>> 238 //// G4cout << " y, xx, yy, zz " << y << " " << xx << " "
>> 239 //// << yy << " " << zz << G4endl;
>> 240 //// G4cout << " phiG, cthetaG, sthetaG " << phiG << " "
>> 241 //// << cthetaG << " "
>> 242 //// << sthetaG << " " << G4endl;
>> 243 //
>> 244 //-----------------------------------------------------------------------
>> 245 //
>> 246 // Maybe certain restrictions on the kinematical variables:
>> 247 //
>> 248 //// if (cthetaE > 0.01)goto leap2;
>> 249 //// if (cthetaG > 0.01)goto leap2;
>> 250 //// if (std::fabs(x-0.5) > 0.5 )goto leap2;
>> 251 //// if (std::fabs(y-0.5) > 0.5 )goto leap2;
>> 252 //
>> 253 //-----------------------------------------------------------------------
>> 254 //
>> 255 // Calculate the angle between positron and photon (cosine)
>> 256 //
>> 257 cthetaGE = cthetaE*cthetaG+sthetaE*sthetaG*std::cos(phiE-phiG);
>> 258 //
>> 259 //// G4cout << x << " " << cthetaE << " " << sthetaE << " "
>> 260 //// << y << " " << cthetaG << " " << sthetaG << " "
>> 261 //// << cthetaGE
>> 262 //
>> 263 //-----------------------------------------------------------------------
>> 264 //
>> 265 G4double term0 = eps*eps;
>> 266 G4double term1 = x*((1.0-eps)*(1.0-eps))+2.0*eps;
>> 267 G4double beta = std::sqrt( x*((1.0-eps)*(1.0-eps))*
>> 268 (x*((1.0-eps)*(1.0-eps))+4.0*eps))/term1;
>> 269 G4double delta = 1.0-beta*cthetaGE;
>> 270
>> 271 G4double term3 = y*(1.0-(eps*eps));
>> 272 G4double term6 = term1*delta*term3;
>> 273
>> 274 G4double Qsqr = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps));
>> 275 //
>> 276 //-----------------------------------------------------------------------
>> 277 //
>> 278 // Check the kinematics.
>> 279 //
>> 280 if ( Qsqr>=0.0 && Qsqr<=1.0 ) break;
>> 281 }
>> 282 //
>> 283 //// G4cout << x << " " << y << " " << beta << " " << Qsqr << G4endl;
>> 284 //
>> 285 // Do the calculation for -1 muon polarization (i.e. mu+)
>> 286 //
>> 287 G4double Pmu = -1.0;
>> 288 if (GetParentName() == "mu-")Pmu = +1.0;
>> 289 //
>> 290 // and for Fronsdal
>> 291 //
>> 292 //-----------------------------------------------------------------------
>> 293 //
>> 294 som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE);
>> 295 //
>> 296 //// if(som0<0.0){
>> 297 //// G4cout << " som0 < 0 in Fronsdal " << som0
>> 298 //// << " at event " << i << G4endl;
>> 299 //// G4cout << Pmu << " " << x << " " << y << " "
>> 300 //// << cthetaE << " " << cthetaG << " "
>> 301 //// << cthetaGE << " " << som0 << G4endl;
>> 302 //// }
>> 303 //
>> 304 //-----------------------------------------------------------------------
>> 305 //
>> 306 //// G4cout << x << " " << y << " " << som0 << G4endl;
>> 307 //
>> 308 //----------------------------------------------------------------------
>> 309 //
>> 310 // Sample the decay rate
>> 311 //
>> 312
>> 313 if (G4UniformRand()*177.0 <= som0) break;
249 } 314 }
250 315
251 G4double E = EMMU / 2. * (x * ((1. - eps) * << 316 /// if(i<10000000)goto leap1:
252 G4double G = EMMU / 2. * y * (1. - eps * eps << 317 //
>> 318 //-----------------------------------------------------------------------
>> 319 //
>> 320 G4double E = EMMU/2.*(x*((1.-eps)*(1.-eps))+2.*eps);
>> 321 G4double G = EMMU/2.*y*(1.-eps*eps);
>> 322 //
>> 323 //-----------------------------------------------------------------------
>> 324 //
253 325
254 if (E < EMASS) E = EMASS; << 326 if(E < EMASS) E = EMASS;
255 327
256 // calculate daughter momentum 328 // calculate daughter momentum
257 G4double daughtermomentum[4]; 329 G4double daughtermomentum[4];
258 330
259 daughtermomentum[0] = std::sqrt(E * E - EMAS << 331 daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS);
260 332
261 G4double sthetaE = std::sqrt(1. - cthetaE * << 333 G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE);
262 G4double cphiE = std::cos(phiE); 334 G4double cphiE = std::cos(phiE);
263 G4double sphiE = std::sin(phiE); 335 G4double sphiE = std::sin(phiE);
264 336
265 // Coordinates of the decay positron with re << 337 //Coordinates of the decay positron with respect to the muon spin
266 338
267 G4double px = sthetaE * cphiE; << 339 G4double px = sthetaE*cphiE;
268 G4double py = sthetaE * sphiE; << 340 G4double py = sthetaE*sphiE;
269 G4double pz = cthetaE; 341 G4double pz = cthetaE;
270 342
271 G4ThreeVector direction0(px, py, pz); << 343 G4ThreeVector direction0(px,py,pz);
272 344
273 direction0.rotateUz(parent_polarization); 345 direction0.rotateUz(parent_polarization);
274 346
275 auto daughterparticle0 = << 347 G4DynamicParticle * daughterparticle0
276 new G4DynamicParticle(G4MT_daughters[0], d << 348 = new G4DynamicParticle( G4MT_daughters[0], daughtermomentum[0]*direction0);
277 349
278 products->PushProducts(daughterparticle0); 350 products->PushProducts(daughterparticle0);
279 351
280 daughtermomentum[1] = G; 352 daughtermomentum[1] = G;
281 353
282 G4double sthetaG = std::sqrt(1. - cthetaG * << 354 G4double sthetaG = std::sqrt(1.-cthetaG*cthetaG);
283 G4double cphiG = std::cos(phiG); 355 G4double cphiG = std::cos(phiG);
284 G4double sphiG = std::sin(phiG); 356 G4double sphiG = std::sin(phiG);
285 357
286 // Coordinates of the decay gamma with respe << 358 //Coordinates of the decay gamma with respect to the muon spin
287 359
288 px = sthetaG * cphiG; << 360 px = sthetaG*cphiG;
289 py = sthetaG * sphiG; << 361 py = sthetaG*sphiG;
290 pz = cthetaG; 362 pz = cthetaG;
291 363
292 G4ThreeVector direction1(px, py, pz); << 364 G4ThreeVector direction1(px,py,pz);
293 365
294 direction1.rotateUz(parent_polarization); 366 direction1.rotateUz(parent_polarization);
295 367
296 auto daughterparticle1 = << 368 G4DynamicParticle * daughterparticle1
297 new G4DynamicParticle(G4MT_daughters[1], d << 369 = new G4DynamicParticle( G4MT_daughters[1], daughtermomentum[1]*direction1);
298 370
299 products->PushProducts(daughterparticle1); 371 products->PushProducts(daughterparticle1);
300 372
301 // daughter 3 ,4 (neutrinos) 373 // daughter 3 ,4 (neutrinos)
302 // create neutrinos in the C.M frame of two 374 // create neutrinos in the C.M frame of two neutrinos
303 375
304 G4double energy2 = parentmass - E - G; << 376 G4double energy2 = parentmass*(1.0 - (x+y)/2.0);
305 377
306 G4ThreeVector P34 = -1. * (daughtermomentum[ << 378 G4double vmass = std::sqrt((energy2-
307 G4double vmass2 = energy2 * energy2 - P34.ma << 379 (daughtermomentum[0]+daughtermomentum[1]))*
308 G4double vmass = std::sqrt(vmass2); << 380 (energy2+
309 << 381 (daughtermomentum[0]+daughtermomentum[1])));
310 G4double costhetan = 2. * G4UniformRand() - << 382 G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2;
311 G4double sinthetan = std::sqrt((1.0 - costhe << 383 beta = -1.0 * std::min(beta,0.99);
312 G4double phin = twopi * G4UniformRand() * ra << 384
>> 385 G4double costhetan = 2.*G4UniformRand()-1.0;
>> 386 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
>> 387 G4double phin = twopi*G4UniformRand()*rad;
313 G4double sinphin = std::sin(phin); 388 G4double sinphin = std::sin(phin);
314 G4double cosphin = std::cos(phin); 389 G4double cosphin = std::cos(phin);
315 390
316 G4ThreeVector direction2(sinthetan * cosphin << 391 G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan);
317 392
318 auto daughterparticle2 = new G4DynamicPartic << 393 G4DynamicParticle * daughterparticle2
319 auto daughterparticle3 = << 394 = new G4DynamicParticle( G4MT_daughters[2], direction2*(vmass/2.));
320 new G4DynamicParticle(G4MT_daughters[3], d << 395 G4DynamicParticle * daughterparticle3
>> 396 = new G4DynamicParticle( G4MT_daughters[3], direction2*(-1.0*vmass/2.));
321 397
322 // boost to the muon rest frame 398 // boost to the muon rest frame
323 G4double beta = P34.mag() / energy2; << 399
324 G4ThreeVector direction34 = P34.unit(); << 400 G4ThreeVector direction34(direction0.x()+direction1.x(),
>> 401 direction0.y()+direction1.y(),
>> 402 direction0.z()+direction1.z());
>> 403 direction34 = direction34.unit();
325 404
326 G4LorentzVector p4 = daughterparticle2->Get4 405 G4LorentzVector p4 = daughterparticle2->Get4Momentum();
327 p4.boost(direction34.x() * beta, direction34 << 406 p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
328 daughterparticle2->Set4Momentum(p4); 407 daughterparticle2->Set4Momentum(p4);
329 408
330 p4 = daughterparticle3->Get4Momentum(); 409 p4 = daughterparticle3->Get4Momentum();
331 p4.boost(direction34.x() * beta, direction34 << 410 p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta);
332 daughterparticle3->Set4Momentum(p4); 411 daughterparticle3->Set4Momentum(p4);
333 412
334 products->PushProducts(daughterparticle2); 413 products->PushProducts(daughterparticle2);
335 products->PushProducts(daughterparticle3); 414 products->PushProducts(daughterparticle3);
336 415
337 daughtermomentum[2] = daughterparticle2->Get 416 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
338 daughtermomentum[3] = daughterparticle3->Get 417 daughtermomentum[3] = daughterparticle3->GetTotalMomentum();
339 418
340 // output message << 419 // output message
341 #ifdef G4VERBOSE 420 #ifdef G4VERBOSE
342 if (GetVerboseLevel() > 1) { << 421 if (GetVerboseLevel()>1) {
343 G4cout << "G4MuonRadiativeDecayChannelWith << 422 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt ";
344 G4cout << " create decay products in rest << 423 G4cout << " create decay products in rest frame " <<G4endl;
345 G4double TT = daughterparticle0->GetTotalE << 424 products->DumpInfo();
346 + daughterparticle2->GetTota <<
347 G4cout << "e :" << daughterparticle0->G <<
348 G4cout << "gamma:" << daughterparticle1->G <<
349 G4cout << "nu2 :" << daughterparticle2->G <<
350 G4cout << "nu2 :" << daughterparticle3->G <<
351 G4cout << "total:" << (TT - parentmass) / <<
352 if (GetVerboseLevel() > 1) { <<
353 products->DumpInfo(); <<
354 } <<
355 } 425 }
356 #endif 426 #endif
357 <<
358 return products; 427 return products;
359 } 428 }
360 429
361 G4double G4MuonRadiativeDecayChannelWithSpin:: << 430 G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu,
362 << 431 G4double x,
>> 432 G4double y,
>> 433 G4double cthetaE,
>> 434 G4double cthetaG,
363 435 G4double cthetaGE)
364 { 436 {
365 G4double mu = 105.65; << 437 G4double mu = 105.65;
366 G4double me = 0.511; << 438 G4double me = 0.511;
367 G4double rho = 0.75; << 439 G4double rho = 0.75;
368 G4double del = 0.75; << 440 G4double del = 0.75;
369 G4double eps = 0.0; << 441 G4double eps = 0.0;
370 G4double kap = 0.0; << 442 G4double kap = 0.0;
371 G4double ksi = 1.0; << 443 G4double ksi = 1.0;
372 << 444
373 G4double delta = 1 - cthetaGE; << 445 G4double delta = 1-cthetaGE;
374 << 446
375 // Calculation of the functions f(x,y) << 447 // Calculation of the functions f(x,y)
376 << 448
377 G4double f_1s = 12.0 << 449 G4double f_1s = 12.0*((y*y)*(1.0-y)+x*y*(2.0-3.0*y)
378 * ((y * y) * (1.0 - y) + x * << 450 +2.0*(x*x)*(1.0-2.0*y)-2.0*(x*x*x));
379 - 2.0 * (x * x * x)); << 451 G4double f0s = 6.0*(-x*y*(2.0-3.0*(y*y))
380 G4double f0s = 6.0 << 452 -2.0*(x*x)*(1.0-y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
381 * (-x * y * (2.0 - 3.0 * (y * << 453 G4double f1s = 3.0*((x*x)*y*(2.0-3.0*y-3.0*(y*y))
382 + 2.0 * (x * x * x) * (1.0 << 454 -(x*x*x)*y*(4.0+3.0*y));
383 G4double f1s = << 455 G4double f2s = 1.5*((x*x*x)*(y*y)*(2.0+y));
384 3.0 * ((x * x) * y * (2.0 - 3.0 * y - 3.0 << 456
385 G4double f2s = 1.5 * ((x * x * x) * (y * y) << 457 G4double f_1se = 12.0*(x*y*(1.0-y)+(x*x)*(2.0-3.0*y)
386 << 458 -2.0*(x*x*x));
387 G4double f_1se = 12.0 * (x * y * (1.0 - y) + << 459 G4double f0se = 6.0*(-(x*x)*(2.0-y-2.0*(y*y))
388 G4double f0se = 6.0 * (-(x * x) * (2.0 - y - << 460 +(x*x*x)*(2.0+3.0*y));
389 G4double f1se = -3.0 * (x * x * x) * y * (2. << 461 G4double f1se = -3.0*(x*x*x)*y*(2.0+y);
390 G4double f2se = 0.0; << 462 G4double f2se = 0.0;
391 << 463
392 G4double f_1sg = 12.0 * ((y * y) * (1.0 - y) << 464 G4double f_1sg = 12.0*((y*y)*(1.0-y)+x*y*(1.0-2.0*y)
393 G4double f0sg = << 465 -(x*x)*y);
394 6.0 * (-x * (y * y) * (2.0 - 3.0 * y) - (x << 466 G4double f0sg = 6.0*(-x*(y*y)*(2.0-3.0*y)-(x*x)*y*(1.0-4.0*y)
395 G4double f1sg = 3.0 * ((x * x) * (y * y) * ( << 467 +(x*x*x)*y);
396 G4double f2sg = 1.5 * (x * x * x) * (y * y * << 468 G4double f1sg = 3.0*((x*x)*(y*y)*(1.0-3.0*y)
397 << 469 -2.0*(x*x*x)*(y*y));
398 G4double f_1v = 8.0 << 470 G4double f2sg = 1.5*(x*x*x)*(y*y*y);
399 * ((y * y) * (3.0 - 2.0 * y) << 471
400 + 2.0 * (x * x) * (3.0 - << 472 G4double f_1v = 8.0*((y*y)*(3.0-2.0*y)+6.0*x*y*(1.0-y)
401 G4double f0v = 8.0 << 473 +2.0*(x*x)*(3.0-4.0*y)-4.0*(x*x*x));
402 * (-x * y * (3.0 - y - (y * y << 474 G4double f0v = 8.0*(-x*y*(3.0-y-(y*y))-(x*x)*(3.0-y-4.0*(y*y))
403 + 2.0 * (x * x * x) * (1.0 << 475 +2.0*(x*x*x)*(1.0+2.0*y));
404 G4double f1v = << 476 G4double f1v = 2.0*((x*x)*y*(6.0-5.0*y-2.0*(y*y))
405 2.0 * ((x * x) * y * (6.0 - 5.0 * y - 2.0 << 477 -2.0*(x*x*x)*y*(4.0+3.0*y));
406 G4double f2v = 2.0 * (x * x * x) * (y * y) * << 478 G4double f2v = 2.0*(x*x*x)*(y*y)*(2.0+y);
407 << 479
408 G4double f_1ve = << 480 G4double f_1ve = 8.0*(x*y*(1.0-2.0*y)
409 8.0 * (x * y * (1.0 - 2.0 * y) + 2.0 * (x << 481 +2.0*(x*x)*(1.0-3.0*y)-4.0*(x*x*x));
410 G4double f0ve = << 482 G4double f0ve = 4.0*(-(x*x)*(2.0-3.0*y-4.0*(y*y))
411 4.0 * (-(x * x) * (2.0 - 3.0 * y - 4.0 * ( << 483 +2.0*(x*x*x)*(2.0+3.0*y));
412 G4double f1ve = -4.0 * (x * x * x) * y * (2. << 484 G4double f1ve = -4.0*(x*x*x)*y*(2.0+y);
413 G4double f2ve = 0.0; << 485 G4double f2ve = 0.0;
414 << 486
415 G4double f_1vg = 8.0 * ((y * y) * (1.0 - 2.0 << 487 G4double f_1vg = 8.0*((y*y)*(1.0-2.0*y)+x*y*(1.0-4.0*y)
416 G4double f0vg = << 488 -2.0*(x*x)*y);
417 4.0 * (2.0 * x * (y * y) * (1.0 + y) - (x << 489 G4double f0vg = 4.0*(2.0*x*(y*y)*(1.0+y)-(x*x)*y*(1.0-4.0*y)
418 G4double f1vg = 2.0 * ((x * x) * (y * y) * ( << 490 +2.0*(x*x*x)*y);
419 G4double f2vg = 2.0 * (x * x * x) * (y * y * << 491 G4double f1vg = 2.0*((x*x)*(y*y)*(1.0-2.0*y)
420 << 492 -4.0*(x*x*x)*(y*y));
421 G4double f_1t = 8.0 << 493 G4double f2vg = 2.0*(x*x*x)*(y*y*y);
422 * ((y * y) * (3.0 - y) + 3.0 << 494
423 - 2.0 * (x * x * x)); << 495 G4double f_1t = 8.0*((y*y)*(3.0-y)+3.0*x*y*(2.0-y)
424 G4double f0t = 4.0 << 496 +2.0*(x*x)*(3.0-2.0*y)-2.0*(x*x*x));
425 * (-x * y * (6.0 + (y * y)) - << 497 G4double f0t = 4.0*(-x*y*(6.0+(y*y))
426 + 2.0 * (x * x * x) * (1.0 << 498 -2.0*(x*x)*(3.0+y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y));
427 G4double f1t = << 499 G4double f1t = 2.0*((x*x)*y*(6.0-5.0*y+(y*y))
428 2.0 * ((x * x) * y * (6.0 - 5.0 * y + (y * << 500 -(x*x*x)*y*(4.0+3.0*y));
429 G4double f2t = (x * x * x) * (y * y) * (2.0 << 501 G4double f2t = (x*x*x)*(y*y)*(2.0+y);
430 << 502
431 G4double f_1te = -8.0 * (x * y * (1.0 + 3.0 << 503 G4double f_1te = -8.0*(x*y*(1.0+3.0*y)+(x*x)*(2.0+3.0*y)
432 G4double f0te = 4.0 * ((x * x) * (2.0 + 3.0 << 504 +2.0*(x*x*x));
433 G4double f1te = -2.0 * (x * x * x) * y * (2. << 505 G4double f0te = 4.0*((x*x)*(2.0+3.0*y+4.0*(y*y))
434 G4double f2te = 0.0; << 506 +(x*x*x)*(2.0+3.0*y));
435 << 507 G4double f1te = -2.0*(x*x*x)*y*(2.0+y);
436 G4double f_1tg = -8.0 * ((y * y) * (1.0 + y) << 508 G4double f2te = 0.0;
437 G4double f0tg = 4.0 * (x * (y * y) * (2.0 - << 509
438 G4double f1tg = -2.0 * ((x * x) * (y * y) * << 510 G4double f_1tg = -8.0*((y*y)*(1.0+y)+x*y+(x*x)*y);
439 G4double f2tg = (x * x * x) * (y * y * y); << 511 G4double f0tg = 4.0*(x*(y*y)*(2.0-y)+(x*x)*y*(1.0+2.0*y)
440 << 512 +(x*x*x)*y);
441 G4double term = delta + 2.0 * (me * me) / (( << 513 G4double f1tg = -2.0*((x*x)*(y*y)*(1.0-y)+2.0*(x*x*x)*y);
442 term = 1.0 / term; << 514 G4double f2tg = (x*x*x)*(y*y*y);
443 << 515
444 G4double nss = term * f_1s + f0s + delta * f << 516 G4double term = delta+2.0*(me*me)/((mu*mu)*(x*x));
445 G4double nv = term * f_1v + f0v + delta * f1 << 517 term = 1.0/term;
446 G4double nt = term * f_1t + f0t + delta * f1 << 518
447 << 519 G4double nss = term*f_1s+f0s+delta*f1s+(delta*delta)*f2s;
448 G4double nse = term * f_1se + f0se + delta * << 520 G4double nv = term*f_1v+f0v+delta*f1v+(delta*delta)*f2v;
449 G4double nve = term * f_1ve + f0ve + delta * << 521 G4double nt = term*f_1t+f0t+delta*f1t+(delta*delta)*f2t;
450 G4double nte = term * f_1te + f0te + delta * << 522
451 << 523 G4double nse = term*f_1se+f0se+delta*f1se+(delta*delta)*f2se;
452 G4double nsg = term * f_1sg + f0sg + delta * << 524 G4double nve = term*f_1ve+f0ve+delta*f1ve+(delta*delta)*f2ve;
453 G4double nvg = term * f_1vg + f0vg + delta * << 525 G4double nte = term*f_1te+f0te+delta*f1te+(delta*delta)*f2te;
454 G4double ntg = term * f_1tg + f0tg + delta * << 526
455 << 527 G4double nsg = term*f_1sg+f0sg+delta*f1sg+(delta*delta)*f2sg;
456 G4double term1 = nv; << 528 G4double nvg = term*f_1vg+f0vg+delta*f1vg+(delta*delta)*f2vg;
457 G4double term2 = 2.0 * nss + nv - nt; << 529 G4double ntg = term*f_1tg+f0tg+delta*f1tg+(delta*delta)*f2tg;
458 G4double term3 = 2.0 * nss - 2.0 * nv + nt; << 530
459 << 531 G4double term1 = nv;
460 G4double term1e = 1.0 / 3.0 * (1.0 - 4.0 / 3 << 532 G4double term2 = 2.0*nss+nv-nt;
461 G4double term2e = 2.0 * nse + 5.0 * nve - nt << 533 G4double term3 = 2.0*nss-2.0*nv+nt;
462 G4double term3e = 2.0 * nse - 2.0 * nve + nt << 534
463 << 535 G4double term1e = 1.0/3.0*(1.0-4.0/3.0*del);
464 G4double term1g = 1.0 / 3.0 * (1.0 - 4.0 / 3 << 536 G4double term2e = 2.0*nse+5.0*nve-nte;
465 G4double term2g = 2.0 * nsg + 5.0 * nvg - nt << 537 G4double term3e = 2.0*nse-2.0*nve+nte;
466 G4double term3g = 2.0 * nsg - 2.0 * nvg + nt << 538
467 << 539 G4double term1g = 1.0/3.0*(1.0-4.0/3.0*del);
468 G4double som00 = term1 + (1.0 - 4.0 / 3.0 * << 540 G4double term2g = 2.0*nsg+5.0*nvg-ntg;
469 G4double som01 = Pmu * ksi << 541 G4double term3g = 2.0*nsg-2.0*nvg+ntg;
470 * (cthetaE * (nve - term1e << 542
471 + cthetaG * (nvg - term1 << 543 G4double som00 = term1+(1.0-4.0/3.0*rho)*term2+eps*term3;
>> 544 G4double som01 = Pmu*ksi*(cthetaE*(nve-term1e*term2e+kap*term3e)
>> 545 +cthetaG*(nvg-term1g*term2g+kap*term3g));
>> 546
>> 547 G4double som0 = (som00+som01)/y;
>> 548 som0 = fine_structure_const/8./(twopi*twopi*twopi)*som0;
472 549
473 G4double som0 = (som00 + som01) / y; << 550 // G4cout << x << " " << y << " " << som00 << " "
474 som0 = fine_structure_const / 8. / (twopi * << 551 // << som01 << " " << som0 << G4endl;
475 552
476 return som0; << 553 return som0;
477 } 554 }
478 555