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 >> 47 #include "Randomize.hh" 41 #include "G4DecayProducts.hh" 48 #include "G4DecayProducts.hh" 42 #include "G4LorentzVector.hh" 49 #include "G4LorentzVector.hh" 43 #include "G4PhysicalConstants.hh" << 44 #include "G4SystemOfUnits.hh" << 45 #include "Randomize.hh" << 46 50 47 G4MuonRadiativeDecayChannelWithSpin::G4MuonRad << 51 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin() 48 const G4String& theParentName, G4double theB << 52 : G4VDecayChannel(), 49 : G4VDecayChannel("Radiative Muon Decay", 1) << 53 parent_polarization() >> 54 { >> 55 } >> 56 >> 57 G4MuonRadiativeDecayChannelWithSpin:: >> 58 G4MuonRadiativeDecayChannelWithSpin(const G4String& theParentName, >> 59 G4double theBR) >> 60 : G4VDecayChannel("Radiative Muon Decay",1), >> 61 parent_polarization() 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= << 83 { 92 { 84 if (this != &right) { << 93 } >> 94 >> 95 G4MuonRadiativeDecayChannelWithSpin::G4MuonRadiativeDecayChannelWithSpin(const G4MuonRadiativeDecayChannelWithSpin &right): >> 96 G4VDecayChannel(right) >> 97 { >> 98 parent_polarization = right.parent_polarization; >> 99 } >> 100 >> 101 G4MuonRadiativeDecayChannelWithSpin & G4MuonRadiativeDecayChannelWithSpin::operator=(const G4MuonRadiativeDecayChannelWithSpin & right) >> 102 { >> 103 if (this != &right) { 85 kinematics_name = right.kinematics_name; 104 kinematics_name = right.kinematics_name; 86 verboseLevel = right.verboseLevel; 105 verboseLevel = right.verboseLevel; 87 rbranch = right.rbranch; 106 rbranch = right.rbranch; 88 107 89 // copy parent name 108 // copy parent name 90 parent_name = new G4String(*right.parent_n 109 parent_name = new G4String(*right.parent_name); 91 110 92 // clear daughters_name array 111 // clear daughters_name array 93 ClearDaughtersName(); 112 ClearDaughtersName(); 94 113 95 // recreate array 114 // recreate array 96 numberOfDaughters = right.numberOfDaughter 115 numberOfDaughters = right.numberOfDaughters; 97 if (numberOfDaughters > 0) { << 116 if ( numberOfDaughters >0 ) { 98 if (daughters_name != nullptr) ClearDaug << 117 if (daughters_name !=0) ClearDaughtersName(); 99 daughters_name = new G4String*[numberOfD 118 daughters_name = new G4String*[numberOfDaughters]; 100 // copy daughters name << 119 //copy daughters name 101 for (G4int index = 0; index < numberOfDa << 120 for (G4int index=0; index < numberOfDaughters; index++) { 102 daughters_name[index] = new G4String(* << 121 daughters_name[index] = new G4String(*right.daughters_name[index]); 103 } 122 } 104 } 123 } 105 parent_polarization = right.parent_polariz 124 parent_polarization = right.parent_polarization; 106 } 125 } 107 return *this; 126 return *this; 108 } 127 } 109 128 110 G4DecayProducts* G4MuonRadiativeDecayChannelWi << 129 >> 130 G4DecayProducts *G4MuonRadiativeDecayChannelWithSpin::DecayIt(G4double) 111 { 131 { >> 132 112 #ifdef G4VERBOSE 133 #ifdef G4VERBOSE 113 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 134 if (GetVerboseLevel()>1) >> 135 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt "; 114 #endif 136 #endif 115 137 116 CheckAndFillParent(); << 138 if (parent == 0) FillParent(); 117 CheckAndFillDaughters(); << 139 if (daughters == 0) FillDaughters(); 118 140 119 // parent mass 141 // parent mass 120 G4double parentmass = G4MT_parent->GetPDGMas << 142 G4double parentmass = parent->GetPDGMass(); 121 143 122 G4double EMMU = parentmass; 144 G4double EMMU = parentmass; 123 145 124 // daughters'mass << 146 //daughters'mass 125 G4double daughtermass[4]; << 147 G4double daughtermass[4]; 126 // G4double sumofdaughtermass = 0.0; << 148 G4double sumofdaughtermass = 0.0; 127 for (G4int index = 0; index < 4; ++index) { << 149 for (G4int index=0; index<4; index++){ 128 daughtermass[index] = G4MT_daughters[index << 150 daughtermass[index] = daughters[index]->GetPDGMass(); 129 // sumofdaughtermass += daughtermass[index << 151 sumofdaughtermass += daughtermass[index]; 130 } 152 } 131 153 132 G4double EMASS = daughtermass[0]; 154 G4double EMASS = daughtermass[0]; 133 155 134 // create parent G4DynamicParticle at rest << 156 //create parent G4DynamicParticle at rest 135 G4ThreeVector dummy; 157 G4ThreeVector dummy; 136 auto parentparticle = new G4DynamicParticle( << 158 G4DynamicParticle * parentparticle = 137 << 159 new G4DynamicParticle( parent, dummy, 0.0); 138 // create G4Decayproducts << 160 //create G4Decayproducts 139 auto products = new G4DecayProducts(*parentp << 161 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 140 delete parentparticle; 162 delete parentparticle; 141 163 142 G4double eps = EMASS / EMMU; << 164 G4int i = 0; >> 165 >> 166 G4double eps = EMASS/EMMU; 143 167 144 G4double som0, x, y, xx, yy, zz; << 168 G4double som0, Qsqr, x, y, xx, yy, zz; 145 G4double cthetaE, cthetaG, cthetaGE, phiE, p 169 G4double cthetaE, cthetaG, cthetaGE, phiE, phiG; 146 const std::size_t MAX_LOOP = 10000; << 147 170 148 for (std::size_t loop_counter1 = 0; loop_cou << 171 do { 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 do { >> 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 << 197 198 // What you get: << 198 phiE = atan4(xx,yy); 199 // << 199 cthetaE = zz/x; 200 // y = photon energy << 200 G4double sthetaE = std::sqrt((xx*xx)+(yy*yy))/x; 201 // phiG = azimutal angle of photon mo << 201 // 202 // cthetaG = cosine of polar angle of ph << 202 // What you get: 203 // sthetaG = sine of polar angle of phot << 203 // 204 // << 204 // x = positron energy 205 //// G4cout << " y, xx, yy, zz " << y < << 205 // phiE = azimutal angle of positron momentum 206 //// << yy < << 206 // cthetaE = cosine of polar angle of positron momentum 207 //// G4cout << " phiG, cthetaG, sthetaG << 207 // sthetaE = sine of polar angle of positron momentum 208 //// << 208 // 209 //// << 209 //// G4cout << " x, xx, yy, zz " << x << " " << xx << " " 210 << 210 //// << yy << " " << zz << G4endl; 211 // Calculate the angle between posi << 211 //// G4cout << " phiE, cthetaE, sthetaE " << phiE << " " 212 // << 212 //// << cthetaE << " " 213 cthetaGE = cthetaE * cthetaG + sthetaE * << 213 //// << sthetaE << " " << G4endl; 214 << 214 // 215 //// G4cout << x << " " << cthetaE << " << 215 //----------------------------------------------------------------------- 216 //// << y << " " << cthetaG << " << 216 // 217 //// << cthetaGE << 217 // For the photon 218 << 218 // 219 G4double term0 = eps * eps; << 219 y = G4UniformRand(); 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 220 244 som0 = fron(Pmu, x, y, cthetaE, cthetaG, c << 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 } 249 } << 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; 250 270 251 G4double E = EMMU / 2. * (x * ((1. - eps) * << 271 G4double term3 = y*(1.0-(eps*eps)); 252 G4double G = EMMU / 2. * y * (1. - eps * eps << 272 G4double term6 = term1*delta*term3; 253 273 254 if (E < EMASS) E = EMASS; << 274 Qsqr = (1.0-term1-term3+term0+0.5*term6)/((1.0-eps)*(1.0-eps)); >> 275 // >> 276 //----------------------------------------------------------------------- >> 277 // >> 278 // Check the kinematics. >> 279 // >> 280 } while ( Qsqr<0.0 || Qsqr>1.0 ); >> 281 // >> 282 //// G4cout << x << " " << y << " " << beta << " " << Qsqr << G4endl; >> 283 // >> 284 // Do the calculation for -1 muon polarization (i.e. mu+) >> 285 // >> 286 G4double Pmu = -1.0; >> 287 if (GetParentName() == "mu-")Pmu = +1.0; >> 288 // >> 289 // and for Fronsdal >> 290 // >> 291 //----------------------------------------------------------------------- >> 292 // >> 293 som0 = fron(Pmu,x,y,cthetaE,cthetaG,cthetaGE); >> 294 // >> 295 //// if(som0<0.0){ >> 296 //// G4cout << " som0 < 0 in Fronsdal " << som0 >> 297 //// << " at event " << i << G4endl; >> 298 //// G4cout << Pmu << " " << x << " " << y << " " >> 299 //// << cthetaE << " " << cthetaG << " " >> 300 //// << cthetaGE << " " << som0 << G4endl; >> 301 //// } >> 302 // >> 303 //----------------------------------------------------------------------- >> 304 // >> 305 //// G4cout << x << " " << y << " " << som0 << G4endl; >> 306 // >> 307 //---------------------------------------------------------------------- >> 308 // >> 309 // Sample the decay rate >> 310 // >> 311 >> 312 } while (G4UniformRand()*177.0 > som0); >> 313 >> 314 /// if(i<10000000)goto leap1: >> 315 // >> 316 //----------------------------------------------------------------------- >> 317 // >> 318 G4double E = EMMU/2.*(x*((1.-eps)*(1.-eps))+2.*eps); >> 319 G4double G = EMMU/2.*y*(1.-eps*eps); >> 320 // >> 321 //----------------------------------------------------------------------- >> 322 // >> 323 >> 324 if(E < EMASS) E = EMASS; 255 325 256 // calculate daughter momentum 326 // calculate daughter momentum 257 G4double daughtermomentum[4]; 327 G4double daughtermomentum[4]; 258 328 259 daughtermomentum[0] = std::sqrt(E * E - EMAS << 329 daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS); 260 330 261 G4double sthetaE = std::sqrt(1. - cthetaE * << 331 G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE); 262 G4double cphiE = std::cos(phiE); 332 G4double cphiE = std::cos(phiE); 263 G4double sphiE = std::sin(phiE); 333 G4double sphiE = std::sin(phiE); 264 334 265 // Coordinates of the decay positron with re << 335 //Coordinates of the decay positron with respect to the muon spin 266 336 267 G4double px = sthetaE * cphiE; << 337 G4double px = sthetaE*cphiE; 268 G4double py = sthetaE * sphiE; << 338 G4double py = sthetaE*sphiE; 269 G4double pz = cthetaE; 339 G4double pz = cthetaE; 270 340 271 G4ThreeVector direction0(px, py, pz); << 341 G4ThreeVector direction0(px,py,pz); 272 342 273 direction0.rotateUz(parent_polarization); 343 direction0.rotateUz(parent_polarization); 274 344 275 auto daughterparticle0 = << 345 G4DynamicParticle * daughterparticle0 276 new G4DynamicParticle(G4MT_daughters[0], d << 346 = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0); 277 347 278 products->PushProducts(daughterparticle0); 348 products->PushProducts(daughterparticle0); 279 349 280 daughtermomentum[1] = G; 350 daughtermomentum[1] = G; 281 351 282 G4double sthetaG = std::sqrt(1. - cthetaG * << 352 G4double sthetaG = std::sqrt(1.-cthetaG*cthetaG); 283 G4double cphiG = std::cos(phiG); 353 G4double cphiG = std::cos(phiG); 284 G4double sphiG = std::sin(phiG); 354 G4double sphiG = std::sin(phiG); 285 355 286 // Coordinates of the decay gamma with respe << 356 //Coordinates of the decay gamma with respect to the muon spin 287 357 288 px = sthetaG * cphiG; << 358 px = sthetaG*cphiG; 289 py = sthetaG * sphiG; << 359 py = sthetaG*sphiG; 290 pz = cthetaG; 360 pz = cthetaG; 291 361 292 G4ThreeVector direction1(px, py, pz); << 362 G4ThreeVector direction1(px,py,pz); 293 363 294 direction1.rotateUz(parent_polarization); 364 direction1.rotateUz(parent_polarization); 295 365 296 auto daughterparticle1 = << 366 G4DynamicParticle * daughterparticle1 297 new G4DynamicParticle(G4MT_daughters[1], d << 367 = new G4DynamicParticle( daughters[1], daughtermomentum[1]*direction1); 298 368 299 products->PushProducts(daughterparticle1); 369 products->PushProducts(daughterparticle1); 300 370 301 // daughter 3 ,4 (neutrinos) 371 // daughter 3 ,4 (neutrinos) 302 // create neutrinos in the C.M frame of two 372 // create neutrinos in the C.M frame of two neutrinos 303 373 304 G4double energy2 = parentmass - E - G; << 374 G4double energy2 = parentmass*(1.0 - (x+y)/2.0); 305 375 306 G4ThreeVector P34 = -1. * (daughtermomentum[ << 376 G4double vmass = std::sqrt((energy2- 307 G4double vmass2 = energy2 * energy2 - P34.ma << 377 (daughtermomentum[0]+daughtermomentum[1]))* 308 G4double vmass = std::sqrt(vmass2); << 378 (energy2+ 309 << 379 (daughtermomentum[0]+daughtermomentum[1]))); 310 G4double costhetan = 2. * G4UniformRand() - << 380 G4double beta = (daughtermomentum[0]+daughtermomentum[1])/energy2; 311 G4double sinthetan = std::sqrt((1.0 - costhe << 381 beta = -1.0 * std::min(beta,0.99); 312 G4double phin = twopi * G4UniformRand() * ra << 382 >> 383 G4double costhetan = 2.*G4UniformRand()-1.0; >> 384 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); >> 385 G4double phin = twopi*G4UniformRand()*rad; 313 G4double sinphin = std::sin(phin); 386 G4double sinphin = std::sin(phin); 314 G4double cosphin = std::cos(phin); 387 G4double cosphin = std::cos(phin); 315 388 316 G4ThreeVector direction2(sinthetan * cosphin << 389 G4ThreeVector direction2(sinthetan*cosphin,sinthetan*sinphin,costhetan); 317 390 318 auto daughterparticle2 = new G4DynamicPartic << 391 G4DynamicParticle * daughterparticle2 319 auto daughterparticle3 = << 392 = new G4DynamicParticle( daughters[2], direction2*(vmass/2.)); 320 new G4DynamicParticle(G4MT_daughters[3], d << 393 G4DynamicParticle * daughterparticle3 >> 394 = new G4DynamicParticle( daughters[3], direction2*(-1.0*vmass/2.)); 321 395 322 // boost to the muon rest frame 396 // boost to the muon rest frame 323 G4double beta = P34.mag() / energy2; << 397 324 G4ThreeVector direction34 = P34.unit(); << 398 G4ThreeVector direction34(direction0.x()+direction1.x(), >> 399 direction0.y()+direction1.y(), >> 400 direction0.z()+direction1.z()); >> 401 direction34 = direction34.unit(); 325 402 326 G4LorentzVector p4 = daughterparticle2->Get4 403 G4LorentzVector p4 = daughterparticle2->Get4Momentum(); 327 p4.boost(direction34.x() * beta, direction34 << 404 p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta); 328 daughterparticle2->Set4Momentum(p4); 405 daughterparticle2->Set4Momentum(p4); 329 406 330 p4 = daughterparticle3->Get4Momentum(); 407 p4 = daughterparticle3->Get4Momentum(); 331 p4.boost(direction34.x() * beta, direction34 << 408 p4.boost(direction34.x()*beta,direction34.y()*beta,direction34.z()*beta); 332 daughterparticle3->Set4Momentum(p4); 409 daughterparticle3->Set4Momentum(p4); 333 410 334 products->PushProducts(daughterparticle2); 411 products->PushProducts(daughterparticle2); 335 products->PushProducts(daughterparticle3); 412 products->PushProducts(daughterparticle3); 336 413 337 daughtermomentum[2] = daughterparticle2->Get 414 daughtermomentum[2] = daughterparticle2->GetTotalMomentum(); 338 daughtermomentum[3] = daughterparticle3->Get 415 daughtermomentum[3] = daughterparticle3->GetTotalMomentum(); 339 416 340 // output message << 417 // output message 341 #ifdef G4VERBOSE 418 #ifdef G4VERBOSE 342 if (GetVerboseLevel() > 1) { << 419 if (GetVerboseLevel()>1) { 343 G4cout << "G4MuonRadiativeDecayChannelWith << 420 G4cout << "G4MuonRadiativeDecayChannelWithSpin::DecayIt "; 344 G4cout << " create decay products in rest << 421 G4cout << " create decay products in rest frame " <<G4endl; 345 G4double TT = daughterparticle0->GetTotalE << 422 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 } 423 } 356 #endif 424 #endif 357 << 358 return products; 425 return products; 359 } 426 } 360 427 361 G4double G4MuonRadiativeDecayChannelWithSpin:: << 428 G4double G4MuonRadiativeDecayChannelWithSpin::fron(G4double Pmu, 362 << 429 G4double x, >> 430 G4double y, >> 431 G4double cthetaE, >> 432 G4double cthetaG, 363 433 G4double cthetaGE) 364 { 434 { 365 G4double mu = 105.65; << 435 G4double mu = 105.65; 366 G4double me = 0.511; << 436 G4double me = 0.511; 367 G4double rho = 0.75; << 437 G4double rho = 0.75; 368 G4double del = 0.75; << 438 G4double del = 0.75; 369 G4double eps = 0.0; << 439 G4double eps = 0.0; 370 G4double kap = 0.0; << 440 G4double kap = 0.0; 371 G4double ksi = 1.0; << 441 G4double ksi = 1.0; 372 << 442 373 G4double delta = 1 - cthetaGE; << 443 G4double delta = 1-cthetaGE; 374 << 444 375 // Calculation of the functions f(x,y) << 445 // Calculation of the functions f(x,y) 376 << 446 377 G4double f_1s = 12.0 << 447 G4double f_1s = 12.0*((y*y)*(1.0-y)+x*y*(2.0-3.0*y) 378 * ((y * y) * (1.0 - y) + x * << 448 +2.0*(x*x)*(1.0-2.0*y)-2.0*(x*x*x)); 379 - 2.0 * (x * x * x)); << 449 G4double f0s = 6.0*(-x*y*(2.0-3.0*(y*y)) 380 G4double f0s = 6.0 << 450 -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 * << 451 G4double f1s = 3.0*((x*x)*y*(2.0-3.0*y-3.0*(y*y)) 382 + 2.0 * (x * x * x) * (1.0 << 452 -(x*x*x)*y*(4.0+3.0*y)); 383 G4double f1s = << 453 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 << 454 385 G4double f2s = 1.5 * ((x * x * x) * (y * y) << 455 G4double f_1se = 12.0*(x*y*(1.0-y)+(x*x)*(2.0-3.0*y) 386 << 456 -2.0*(x*x*x)); 387 G4double f_1se = 12.0 * (x * y * (1.0 - y) + << 457 G4double f0se = 6.0*(-(x*x)*(2.0-y-2.0*(y*y)) 388 G4double f0se = 6.0 * (-(x * x) * (2.0 - y - << 458 +(x*x*x)*(2.0+3.0*y)); 389 G4double f1se = -3.0 * (x * x * x) * y * (2. << 459 G4double f1se = -3.0*(x*x*x)*y*(2.0+y); 390 G4double f2se = 0.0; << 460 G4double f2se = 0.0; 391 << 461 392 G4double f_1sg = 12.0 * ((y * y) * (1.0 - y) << 462 G4double f_1sg = 12.0*((y*y)*(1.0-y)+x*y*(1.0-2.0*y) 393 G4double f0sg = << 463 -(x*x)*y); 394 6.0 * (-x * (y * y) * (2.0 - 3.0 * y) - (x << 464 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) * ( << 465 +(x*x*x)*y); 396 G4double f2sg = 1.5 * (x * x * x) * (y * y * << 466 G4double f1sg = 3.0*((x*x)*(y*y)*(1.0-3.0*y) 397 << 467 -2.0*(x*x*x)*(y*y)); 398 G4double f_1v = 8.0 << 468 G4double f2sg = 1.5*(x*x*x)*(y*y*y); 399 * ((y * y) * (3.0 - 2.0 * y) << 469 400 + 2.0 * (x * x) * (3.0 - << 470 G4double f_1v = 8.0*((y*y)*(3.0-2.0*y)+6.0*x*y*(1.0-y) 401 G4double f0v = 8.0 << 471 +2.0*(x*x)*(3.0-4.0*y)-4.0*(x*x*x)); 402 * (-x * y * (3.0 - y - (y * y << 472 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 << 473 +2.0*(x*x*x)*(1.0+2.0*y)); 404 G4double f1v = << 474 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 << 475 -2.0*(x*x*x)*y*(4.0+3.0*y)); 406 G4double f2v = 2.0 * (x * x * x) * (y * y) * << 476 G4double f2v = 2.0*(x*x*x)*(y*y)*(2.0+y); 407 << 477 408 G4double f_1ve = << 478 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 << 479 +2.0*(x*x)*(1.0-3.0*y)-4.0*(x*x*x)); 410 G4double f0ve = << 480 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 * ( << 481 +2.0*(x*x*x)*(2.0+3.0*y)); 412 G4double f1ve = -4.0 * (x * x * x) * y * (2. << 482 G4double f1ve = -4.0*(x*x*x)*y*(2.0+y); 413 G4double f2ve = 0.0; << 483 G4double f2ve = 0.0; 414 << 484 415 G4double f_1vg = 8.0 * ((y * y) * (1.0 - 2.0 << 485 G4double f_1vg = 8.0*((y*y)*(1.0-2.0*y)+x*y*(1.0-4.0*y) 416 G4double f0vg = << 486 -2.0*(x*x)*y); 417 4.0 * (2.0 * x * (y * y) * (1.0 + y) - (x << 487 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) * ( << 488 +2.0*(x*x*x)*y); 419 G4double f2vg = 2.0 * (x * x * x) * (y * y * << 489 G4double f1vg = 2.0*((x*x)*(y*y)*(1.0-2.0*y) 420 << 490 -4.0*(x*x*x)*(y*y)); 421 G4double f_1t = 8.0 << 491 G4double f2vg = 2.0*(x*x*x)*(y*y*y); 422 * ((y * y) * (3.0 - y) + 3.0 << 492 423 - 2.0 * (x * x * x)); << 493 G4double f_1t = 8.0*((y*y)*(3.0-y)+3.0*x*y*(2.0-y) 424 G4double f0t = 4.0 << 494 +2.0*(x*x)*(3.0-2.0*y)-2.0*(x*x*x)); 425 * (-x * y * (6.0 + (y * y)) - << 495 G4double f0t = 4.0*(-x*y*(6.0+(y*y)) 426 + 2.0 * (x * x * x) * (1.0 << 496 -2.0*(x*x)*(3.0+y-3.0*(y*y))+2.0*(x*x*x)*(1.0+2.0*y)); 427 G4double f1t = << 497 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 * << 498 -(x*x*x)*y*(4.0+3.0*y)); 429 G4double f2t = (x * x * x) * (y * y) * (2.0 << 499 G4double f2t = (x*x*x)*(y*y)*(2.0+y); 430 << 500 431 G4double f_1te = -8.0 * (x * y * (1.0 + 3.0 << 501 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 << 502 +2.0*(x*x*x)); 433 G4double f1te = -2.0 * (x * x * x) * y * (2. << 503 G4double f0te = 4.0*((x*x)*(2.0+3.0*y+4.0*(y*y)) 434 G4double f2te = 0.0; << 504 +(x*x*x)*(2.0+3.0*y)); 435 << 505 G4double f1te = -2.0*(x*x*x)*y*(2.0+y); 436 G4double f_1tg = -8.0 * ((y * y) * (1.0 + y) << 506 G4double f2te = 0.0; 437 G4double f0tg = 4.0 * (x * (y * y) * (2.0 - << 507 438 G4double f1tg = -2.0 * ((x * x) * (y * y) * << 508 G4double f_1tg = -8.0*((y*y)*(1.0+y)+x*y+(x*x)*y); 439 G4double f2tg = (x * x * x) * (y * y * y); << 509 G4double f0tg = 4.0*(x*(y*y)*(2.0-y)+(x*x)*y*(1.0+2.0*y) 440 << 510 +(x*x*x)*y); 441 G4double term = delta + 2.0 * (me * me) / (( << 511 G4double f1tg = -2.0*((x*x)*(y*y)*(1.0-y)+2.0*(x*x*x)*y); 442 term = 1.0 / term; << 512 G4double f2tg = (x*x*x)*(y*y*y); 443 << 513 444 G4double nss = term * f_1s + f0s + delta * f << 514 G4double term = delta+2.0*(me*me)/((mu*mu)*(x*x)); 445 G4double nv = term * f_1v + f0v + delta * f1 << 515 term = 1.0/term; 446 G4double nt = term * f_1t + f0t + delta * f1 << 516 447 << 517 G4double ns = term*f_1s+f0s+delta*f1s+(delta*delta)*f2s; 448 G4double nse = term * f_1se + f0se + delta * << 518 G4double nv = term*f_1v+f0v+delta*f1v+(delta*delta)*f2v; 449 G4double nve = term * f_1ve + f0ve + delta * << 519 G4double nt = term*f_1t+f0t+delta*f1t+(delta*delta)*f2t; 450 G4double nte = term * f_1te + f0te + delta * << 520 451 << 521 G4double nse = term*f_1se+f0se+delta*f1se+(delta*delta)*f2se; 452 G4double nsg = term * f_1sg + f0sg + delta * << 522 G4double nve = term*f_1ve+f0ve+delta*f1ve+(delta*delta)*f2ve; 453 G4double nvg = term * f_1vg + f0vg + delta * << 523 G4double nte = term*f_1te+f0te+delta*f1te+(delta*delta)*f2te; 454 G4double ntg = term * f_1tg + f0tg + delta * << 524 455 << 525 G4double nsg = term*f_1sg+f0sg+delta*f1sg+(delta*delta)*f2sg; 456 G4double term1 = nv; << 526 G4double nvg = term*f_1vg+f0vg+delta*f1vg+(delta*delta)*f2vg; 457 G4double term2 = 2.0 * nss + nv - nt; << 527 G4double ntg = term*f_1tg+f0tg+delta*f1tg+(delta*delta)*f2tg; 458 G4double term3 = 2.0 * nss - 2.0 * nv + nt; << 528 459 << 529 G4double term1 = nv; 460 G4double term1e = 1.0 / 3.0 * (1.0 - 4.0 / 3 << 530 G4double term2 = 2.0*ns+nv-nt; 461 G4double term2e = 2.0 * nse + 5.0 * nve - nt << 531 G4double term3 = 2.0*ns-2.0*nv+nt; 462 G4double term3e = 2.0 * nse - 2.0 * nve + nt << 532 463 << 533 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 << 534 G4double term2e = 2.0*nse+5.0*nve-nte; 465 G4double term2g = 2.0 * nsg + 5.0 * nvg - nt << 535 G4double term3e = 2.0*nse-2.0*nve+nte; 466 G4double term3g = 2.0 * nsg - 2.0 * nvg + nt << 536 467 << 537 G4double term1g = 1.0/3.0*(1.0-4.0/3.0*del); 468 G4double som00 = term1 + (1.0 - 4.0 / 3.0 * << 538 G4double term2g = 2.0*nsg+5.0*nvg-ntg; 469 G4double som01 = Pmu * ksi << 539 G4double term3g = 2.0*nsg-2.0*nvg+ntg; 470 * (cthetaE * (nve - term1e << 540 471 + cthetaG * (nvg - term1 << 541 G4double som00 = term1+(1.0-4.0/3.0*rho)*term2+eps*term3; >> 542 G4double som01 = Pmu*ksi*(cthetaE*(nve-term1e*term2e+kap*term3e) >> 543 +cthetaG*(nvg-term1g*term2g+kap*term3g)); >> 544 >> 545 G4double som0 = (som00+som01)/y; >> 546 som0 = fine_structure_const/8./(twopi*twopi*twopi)*som0; 472 547 473 G4double som0 = (som00 + som01) / y; << 548 // G4cout << x << " " << y << " " << som00 << " " 474 som0 = fine_structure_const / 8. / (twopi * << 549 // << som01 << " " << som0 << G4endl; 475 550 476 return som0; << 551 return som0; 477 } 552 } 478 553