| 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 // G4MuonDecayChannelWithSpin class implementa 26 // G4MuonDecayChannelWithSpin class implementation
27 // 27 //
28 // References: 28 // References:
29 // - Florian Scheck "Muon Physics", in Physics 29 // - Florian Scheck "Muon Physics", in Physics Reports
30 // (Review Section of Physics Letters) 44, N 30 // (Review Section of Physics Letters) 44, No. 4 (1978)
31 // 187-248. North-Holland Publishing Company 31 // 187-248. North-Holland Publishing Company, Amsterdam at page 210 cc.
32 // - W.E. Fisher and F. Scheck, Nucl. Phys. B8 32 // - W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
33 33
34 // Authors: P.Gumplinger and T.MacPhail, 17 Au << 34 // Authors: P.Gumplinger and T.MacPhail, 17 August 2004
35 // ------------------------------------------- 35 // --------------------------------------------------------------------
36 36
37 #include "G4MuonDecayChannelWithSpin.hh" 37 #include "G4MuonDecayChannelWithSpin.hh"
38 38
39 #include "G4DecayProducts.hh" <<
40 #include "G4LorentzVector.hh" <<
41 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
43 #include "Randomize.hh" 41 #include "Randomize.hh"
44 42
45 G4MuonDecayChannelWithSpin::G4MuonDecayChannel << 43 #include "G4DecayProducts.hh"
46 << 44 #include "G4LorentzVector.hh"
47 : G4MuonDecayChannel(theParentName, theBR) <<
48 {} <<
49 45
50 G4MuonDecayChannelWithSpin& << 46 G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin()
51 G4MuonDecayChannelWithSpin::operator=(const G4 << 47 : G4MuonDecayChannel()
52 { 48 {
53 if (this != &right) { << 49 }
>> 50
>> 51 G4MuonDecayChannelWithSpin::
>> 52 G4MuonDecayChannelWithSpin(const G4String& theParentName,
>> 53 G4double theBR)
>> 54 : G4MuonDecayChannel(theParentName,theBR)
>> 55 {
>> 56 }
>> 57
>> 58 G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin()
>> 59 {
>> 60 }
>> 61
>> 62 G4MuonDecayChannelWithSpin::
>> 63 G4MuonDecayChannelWithSpin(const G4MuonDecayChannelWithSpin& right)
>> 64 : G4MuonDecayChannel(right)
>> 65 {
>> 66 }
>> 67
>> 68 G4MuonDecayChannelWithSpin& G4MuonDecayChannelWithSpin::
>> 69 operator=(const G4MuonDecayChannelWithSpin& right)
>> 70 {
>> 71 if (this != &right)
>> 72 {
54 kinematics_name = right.kinematics_name; 73 kinematics_name = right.kinematics_name;
55 verboseLevel = right.verboseLevel; 74 verboseLevel = right.verboseLevel;
56 rbranch = right.rbranch; 75 rbranch = right.rbranch;
57 76
58 // copy parent name 77 // copy parent name
59 delete parent_name; 78 delete parent_name;
60 parent_name = new G4String(*right.parent_n 79 parent_name = new G4String(*right.parent_name);
61 80
62 // clear daughters_name array 81 // clear daughters_name array
63 ClearDaughtersName(); 82 ClearDaughtersName();
64 83
65 // recreate array 84 // recreate array
66 numberOfDaughters = right.numberOfDaughter 85 numberOfDaughters = right.numberOfDaughters;
67 if (numberOfDaughters > 0) { << 86 if ( numberOfDaughters > 0 )
>> 87 {
68 daughters_name = new G4String*[numberOfD 88 daughters_name = new G4String*[numberOfDaughters];
69 // copy daughters name 89 // copy daughters name
70 for (G4int index = 0; index < numberOfDa << 90 for (G4int index=0; index<numberOfDaughters; ++index)
>> 91 {
71 daughters_name[index] = new G4String(* 92 daughters_name[index] = new G4String(*right.daughters_name[index]);
72 } 93 }
73 } 94 }
74 } 95 }
75 return *this; 96 return *this;
76 } 97 }
77 98
78 G4DecayProducts* G4MuonDecayChannelWithSpin::D << 99 G4DecayProducts* G4MuonDecayChannelWithSpin::DecayIt(G4double)
79 { 100 {
80 // This version assumes V-A coupling with 1s 101 // This version assumes V-A coupling with 1st order radiative correctons,
81 // the standard model Michel pa 102 // the standard model Michel parameter values, but
82 // gives incorrect energy spect 103 // gives incorrect energy spectrum for neutrinos
83 104
84 #ifdef G4VERBOSE 105 #ifdef G4VERBOSE
85 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 106 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
86 #endif 107 #endif
87 108
88 CheckAndFillParent(); 109 CheckAndFillParent();
89 CheckAndFillDaughters(); 110 CheckAndFillDaughters();
90 111
91 // parent mass 112 // parent mass
92 G4double parentmass = G4MT_parent->GetPDGMas 113 G4double parentmass = G4MT_parent->GetPDGMass();
93 114
94 G4double EMMU = parentmass; 115 G4double EMMU = parentmass;
95 116
96 // daughters'mass << 117 //daughters'mass
97 G4double daughtermass[3]; << 118 G4double daughtermass[3];
98 // G4double sumofdaughtermass = 0.0; << 119 G4double sumofdaughtermass = 0.0;
99 for (G4int index = 0; index < 3; ++index) { << 120 for (G4int index=0; index<3; ++index)
>> 121 {
100 daughtermass[index] = G4MT_daughters[index 122 daughtermass[index] = G4MT_daughters[index]->GetPDGMass();
101 // sumofdaughtermass += daughtermass[index << 123 sumofdaughtermass += daughtermass[index];
102 } 124 }
103 125
104 G4double EMASS = daughtermass[0]; 126 G4double EMASS = daughtermass[0];
105 127
106 // create parent G4DynamicParticle at rest 128 // create parent G4DynamicParticle at rest
107 G4ThreeVector dummy; 129 G4ThreeVector dummy;
108 auto parentparticle = new G4DynamicParticle( << 130 G4DynamicParticle* parentparticle
>> 131 = new G4DynamicParticle( G4MT_parent, dummy, 0.0);
109 // create G4Decayproducts 132 // create G4Decayproducts
110 auto products = new G4DecayProducts(*parentp << 133 G4DecayProducts *products = new G4DecayProducts(*parentparticle);
111 delete parentparticle; 134 delete parentparticle;
112 135
113 // calculate electron energy 136 // calculate electron energy
114 137
115 G4double michel_rho = 0.75; // Standard Mod << 138 G4double michel_rho = 0.75; //Standard Model Michel rho
116 G4double michel_delta = 0.75; // Standard M << 139 G4double michel_delta = 0.75; //Standard Model Michel delta
117 G4double michel_xsi = 1.00; // Standard Mod << 140 G4double michel_xsi = 1.00; //Standard Model Michel xsi
118 G4double michel_eta = 0.00; // Standard Mod << 141 G4double michel_eta = 0.00; //Standard Model eta
119 142
120 G4double rndm, x, ctheta; 143 G4double rndm, x, ctheta;
121 144
122 G4double FG; << 145 G4double FG;
123 G4double FG_max = 2.00; 146 G4double FG_max = 2.00;
124 147
125 G4double W_mue = (EMMU * EMMU + EMASS * EMAS << 148 G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
126 G4double x0 = EMASS / W_mue; << 149 G4double x0 = EMASS/W_mue;
127 <<
128 G4double x0_squared = x0 * x0; <<
129 150
>> 151 G4double x0_squared = x0*x0;
>> 152
130 // ***************************************** 153 // ***************************************************
131 // x0 <= x <= 1. and -1 <= y <= 1 154 // x0 <= x <= 1. and -1 <= y <= 1
132 // 155 //
133 // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g 156 // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y
134 // ***************************************** 157 // ***************************************************
135 158
136 // ***** sampling F(x,y) directly (brute for 159 // ***** sampling F(x,y) directly (brute force) *****
137 160
138 const std::size_t MAX_LOOP = 10000; << 161 const std::size_t MAX_LOOP=10000;
139 for (std::size_t loop_count = 0; loop_count << 162 for (std::size_t loop_count=0; loop_count<MAX_LOOP; ++loop_count)
>> 163 {
140 // Sample the positron energy by sampling 164 // Sample the positron energy by sampling from F
141 165
142 rndm = G4UniformRand(); 166 rndm = G4UniformRand();
143 167
144 x = x0 + rndm * (1. - x0); << 168 x = x0 + rndm*(1.-x0);
145 169
146 G4double x_squared = x * x; << 170 G4double x_squared = x*x;
147 171
148 G4double F_IS, F_AS, G_IS, G_AS; 172 G4double F_IS, F_AS, G_IS, G_AS;
149 173
150 F_IS = 1. / 6. * (-2. * x_squared + 3. * x << 174 F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
151 F_AS = 1. / 6. * std::sqrt(x_squared - x0_ << 175 F_AS = 1./6.*std::sqrt(x_squared-x0_squared)
152 << 176 *(2.*x-2.+std::sqrt(1.-x0_squared));
153 G_IS = 2. / 9. * (michel_rho - 0.75) * (4. << 177
154 G_IS = G_IS + michel_eta * (1. - x) * x0; << 178 G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
155 << 179 G_IS = G_IS + michel_eta*(1.-x)*x0;
156 G_AS = 3. * (michel_xsi - 1.) * (1. - x); << 180
157 G_AS = << 181 G_AS = 3.*(michel_xsi-1.)*(1.-x);
158 G_AS + 2. * (michel_xsi * michel_delta - << 182 G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)
159 G_AS = 1. / 9. * std::sqrt(x_squared - x0_ << 183 *(4.*x-4.+std::sqrt(1.-x0_squared));
>> 184 G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
160 185
161 F_IS = F_IS + G_IS; 186 F_IS = F_IS + G_IS;
162 F_AS = F_AS + G_AS; 187 F_AS = F_AS + G_AS;
163 188
164 // *** Radiative Corrections *** 189 // *** Radiative Corrections ***
165 190
166 const G4double omega = std::log(EMMU / EMA << 191 const G4double omega = std::log(EMMU/EMASS);
167 G4double R_IS = F_c(x, x0, omega); << 192 G4double R_IS = F_c(x,x0,omega);
168 193
169 G4double F = 6. * F_IS + R_IS / std::sqrt( << 194 G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
170 195
171 // *** Radiative Corrections *** 196 // *** Radiative Corrections ***
172 197
173 G4double R_AS = F_theta(x, x0, omega); << 198 G4double R_AS = F_theta(x,x0,omega);
174 199
175 rndm = G4UniformRand(); 200 rndm = G4UniformRand();
176 201
177 ctheta = 2. * rndm - 1.; << 202 ctheta = 2.*rndm-1.;
178 203
179 G4double G = 6. * F_AS - R_AS / std::sqrt( << 204 G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
180 205
181 FG = std::sqrt(x_squared - x0_squared) * F << 206 FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
182 207
183 if (FG > FG_max) { << 208 if(FG>FG_max)
184 G4Exception("G4MuonDecayChannelWithSpin: << 209 {
>> 210 G4Exception("G4MuonDecayChannelWithSpin::DecayIt()",
>> 211 "PART113", JustWarning,
185 "Problem in Muon Decay: FG > 212 "Problem in Muon Decay: FG > FG_max");
186 FG_max = FG; 213 FG_max = FG;
187 } 214 }
188 215
189 rndm = G4UniformRand(); 216 rndm = G4UniformRand();
190 217
191 if (FG >= rndm * FG_max) break; << 218 if (FG >= rndm*FG_max) break;
192 } 219 }
193 220
194 G4double energy = x * W_mue; 221 G4double energy = x * W_mue;
195 222
196 rndm = G4UniformRand(); 223 rndm = G4UniformRand();
197 224
198 G4double phi = twopi * rndm; 225 G4double phi = twopi * rndm;
199 226
200 if (energy < EMASS) energy = EMASS; << 227 if(energy < EMASS) energy = EMASS;
201 228
202 // Calculate daughter momentum 229 // Calculate daughter momentum
203 G4double daughtermomentum[3]; 230 G4double daughtermomentum[3];
204 231
205 daughtermomentum[0] = std::sqrt(energy * ene << 232 daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
206 233
207 G4double stheta = std::sqrt(1. - ctheta * ct << 234 G4double stheta = std::sqrt(1.-ctheta*ctheta);
208 G4double cphi = std::cos(phi); 235 G4double cphi = std::cos(phi);
209 G4double sphi = std::sin(phi); 236 G4double sphi = std::sin(phi);
210 237
211 // Coordinates of the decay positron with re 238 // Coordinates of the decay positron with respect to the muon spin
212 G4double px = stheta * cphi; << 239 G4double px = stheta*cphi;
213 G4double py = stheta * sphi; << 240 G4double py = stheta*sphi;
214 G4double pz = ctheta; 241 G4double pz = ctheta;
215 242
216 G4ThreeVector direction0(px, py, pz); << 243 G4ThreeVector direction0(px,py,pz);
217 244
218 direction0.rotateUz(parent_polarization); 245 direction0.rotateUz(parent_polarization);
219 246
220 auto daughterparticle0 = << 247 G4DynamicParticle * daughterparticle0
221 new G4DynamicParticle(G4MT_daughters[0], d << 248 = new G4DynamicParticle( G4MT_daughters[0], daughtermomentum[0]*direction0);
222 249
223 products->PushProducts(daughterparticle0); 250 products->PushProducts(daughterparticle0);
224 251
>> 252
225 // daughter 1 ,2 (neutrinos) 253 // daughter 1 ,2 (neutrinos)
226 // create neutrinos in the C.M frame of two 254 // create neutrinos in the C.M frame of two neutrinos
227 G4double energy2 = parentmass - energy; << 255 G4double energy2 = parentmass-energy;
228 G4double vmass = std::sqrt((energy2 - daught << 256 G4double vmass = std::sqrt((energy2-daughtermomentum[0])
229 G4double beta = -1.0 * daughtermomentum[0] / << 257 * (energy2+daughtermomentum[0]));
230 G4double costhetan = 2. * G4UniformRand() - << 258 G4double beta = -1.0*daughtermomentum[0]/energy2;
231 G4double sinthetan = std::sqrt((1.0 - costhe << 259 G4double costhetan = 2.*G4UniformRand()-1.0;
232 G4double phin = twopi * G4UniformRand() * ra << 260 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
>> 261 G4double phin = twopi*G4UniformRand()*rad;
233 G4double sinphin = std::sin(phin); 262 G4double sinphin = std::sin(phin);
234 G4double cosphin = std::cos(phin); 263 G4double cosphin = std::cos(phin);
235 264
236 G4ThreeVector direction1(sinthetan * cosphin << 265 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
237 auto daughterparticle1 = new G4DynamicPartic << 266 G4DynamicParticle * daughterparticle1
238 auto daughterparticle2 = << 267 = new G4DynamicParticle( G4MT_daughters[1], direction1*(vmass/2.));
239 new G4DynamicParticle(G4MT_daughters[2], d << 268 G4DynamicParticle * daughterparticle2
>> 269 = new G4DynamicParticle( G4MT_daughters[2], direction1*(-1.0*vmass/2.));
240 270
241 // boost to the muon rest frame 271 // boost to the muon rest frame
242 G4LorentzVector p4; 272 G4LorentzVector p4;
243 p4 = daughterparticle1->Get4Momentum(); 273 p4 = daughterparticle1->Get4Momentum();
244 p4.boost(direction0.x() * beta, direction0.y << 274 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
245 daughterparticle1->Set4Momentum(p4); 275 daughterparticle1->Set4Momentum(p4);
246 p4 = daughterparticle2->Get4Momentum(); 276 p4 = daughterparticle2->Get4Momentum();
247 p4.boost(direction0.x() * beta, direction0.y << 277 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
248 daughterparticle2->Set4Momentum(p4); 278 daughterparticle2->Set4Momentum(p4);
249 products->PushProducts(daughterparticle1); 279 products->PushProducts(daughterparticle1);
250 products->PushProducts(daughterparticle2); 280 products->PushProducts(daughterparticle2);
251 daughtermomentum[1] = daughterparticle1->Get 281 daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
252 daughtermomentum[2] = daughterparticle2->Get 282 daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
253 283
254 // output message 284 // output message
255 #ifdef G4VERBOSE 285 #ifdef G4VERBOSE
256 if (GetVerboseLevel() > 1) { << 286 if (GetVerboseLevel()>1)
>> 287 {
257 G4cout << "G4MuonDecayChannelWithSpin::Dec 288 G4cout << "G4MuonDecayChannelWithSpin::DecayIt ";
258 G4cout << " create decay products in rest << 289 G4cout << " create decay products in rest frame " <<G4endl;
259 G4double TT = daughterparticle0->GetTotalE << 290 G4double TT = daughterparticle0->GetTotalEnergy()
260 + daughterparticle2->GetTota << 291 + daughterparticle1->GetTotalEnergy()
261 G4cout << "e " << daughterparticle0->GetT << 292 + daughterparticle2->GetTotalEnergy();
262 G4cout << "nu1" << daughterparticle1->GetT << 293 G4cout << "e " << daughterparticle0->GetTotalEnergy()/MeV << G4endl;
263 G4cout << "nu2" << daughterparticle2->GetT << 294 G4cout << "nu1" << daughterparticle1->GetTotalEnergy()/MeV << G4endl;
264 G4cout << "total" << (TT - parentmass) / k << 295 G4cout << "nu2" << daughterparticle2->GetTotalEnergy()/MeV << G4endl;
265 if (GetVerboseLevel() > 2) { << 296 G4cout << "total" << (TT-parentmass)/keV << G4endl;
266 products->DumpInfo(); << 297 if (GetVerboseLevel()>2) { products->DumpInfo(); }
267 } <<
268 } 298 }
269 #endif 299 #endif
270 300
271 return products; 301 return products;
272 } 302 }
273 303
274 G4double G4MuonDecayChannelWithSpin::R_c(G4dou << 304 G4double G4MuonDecayChannelWithSpin::R_c(G4double x,G4double omega)
275 { 305 {
276 auto n_max = (G4int)(100. * x); << 306 G4int n_max = (G4int)(100.*x);
277 307
278 if (n_max < 10) n_max = 10; << 308 if(n_max<10)n_max=10;
279 309
280 G4double L2 = 0.0; 310 G4double L2 = 0.0;
281 311
282 for (G4int n = 1; n <= n_max; ++n) { << 312 for(G4int n=1; n<=n_max; ++n)
283 L2 += std::pow(x, n) / (n * n); << 313 {
>> 314 L2 += std::pow(x,n)/(n*n);
284 } 315 }
285 316
286 G4double r_c; 317 G4double r_c;
287 318
288 r_c = 2. * L2 - (pi * pi / 3.) - 2.; << 319 r_c = 2.*L2-(pi*pi/3.)-2.;
289 r_c = r_c + omega * (1.5 + 2. * std::log((1. << 320 r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x));
290 r_c = r_c - std::log(x) * (2. * std::log(x) << 321 r_c = r_c - std::log(x)*(2.*std::log(x)-1.);
291 r_c = r_c + (3. * std::log(x) - 1. - 1. / x) << 322 r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x);
292 323
293 return r_c; 324 return r_c;
294 } 325 }
295 326