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