Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4MuonDecayChannelWithSpin class implementa << 23 // ------------------------------------------------------------ >> 24 // GEANT 4 class header file >> 25 // >> 26 // History: >> 27 // 17 August 2004 P.Gumplinger and T.MacPhail >> 28 // samples Michel spectrum including 1st order >> 29 // radiative corrections >> 30 // Reference: Florian Scheck "Muon Physics", in Physics Reports >> 31 // (Review Section of Physics Letters) 44, No. 4 (1978) >> 32 // 187-248. North-Holland Publishing Company, Amsterdam >> 33 // at page 210 cc. >> 34 // >> 35 // W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25. >> 36 // >> 37 // ------------------------------------------------------------ 27 // 38 // 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" 39 #include "G4MuonDecayChannelWithSpin.hh" 38 40 >> 41 #include "Randomize.hh" 39 #include "G4DecayProducts.hh" 42 #include "G4DecayProducts.hh" 40 #include "G4LorentzVector.hh" 43 #include "G4LorentzVector.hh" 41 #include "G4PhysicalConstants.hh" << 42 #include "G4SystemOfUnits.hh" << 43 #include "Randomize.hh" << 44 44 45 G4MuonDecayChannelWithSpin::G4MuonDecayChannel << 45 G4MuonDecayChannelWithSpin::G4MuonDecayChannelWithSpin(const G4String& theParentName, 46 << 46 G4double theBR) 47 : G4MuonDecayChannel(theParentName, theBR) << 47 : G4MuonDecayChannel(theParentName,theBR) 48 {} << 48 { >> 49 } 49 50 50 G4MuonDecayChannelWithSpin& << 51 G4MuonDecayChannelWithSpin::~G4MuonDecayChannelWithSpin() 51 G4MuonDecayChannelWithSpin::operator=(const G4 << 52 { 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 } 53 } 77 54 78 G4DecayProducts* G4MuonDecayChannelWithSpin::D << 55 G4DecayProducts *G4MuonDecayChannelWithSpin::DecayIt(G4double) 79 { 56 { 80 // This version assumes V-A coupling with 1s 57 // This version assumes V-A coupling with 1st order radiative correctons, 81 // the standard model Michel pa 58 // the standard model Michel parameter values, but 82 // gives incorrect energy spect 59 // gives incorrect energy spectrum for neutrinos 83 60 84 #ifdef G4VERBOSE 61 #ifdef G4VERBOSE 85 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 62 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannelWithSpin::DecayIt "; 86 #endif 63 #endif 87 64 88 CheckAndFillParent(); << 65 if (parent == 0) FillParent(); 89 CheckAndFillDaughters(); << 66 if (daughters == 0) FillDaughters(); 90 67 91 // parent mass 68 // parent mass 92 G4double parentmass = G4MT_parent->GetPDGMas << 69 G4double parentmass = parent->GetPDGMass(); 93 70 94 G4double EMMU = parentmass; << 71 EMMU = parentmass; 95 72 96 // daughters'mass << 73 //daughters'mass 97 G4double daughtermass[3]; << 74 G4double daughtermass[3]; 98 // G4double sumofdaughtermass = 0.0; << 75 G4double sumofdaughtermass = 0.0; 99 for (G4int index = 0; index < 3; ++index) { << 76 for (G4int index=0; index<3; index++){ 100 daughtermass[index] = G4MT_daughters[index << 77 daughtermass[index] = daughters[index]->GetPDGMass(); 101 // sumofdaughtermass += daughtermass[index << 78 sumofdaughtermass += daughtermass[index]; 102 } 79 } 103 80 104 G4double EMASS = daughtermass[0]; << 81 EMASS = daughtermass[0]; 105 82 106 // create parent G4DynamicParticle at rest << 83 //create parent G4DynamicParticle at rest 107 G4ThreeVector dummy; 84 G4ThreeVector dummy; 108 auto parentparticle = new G4DynamicParticle( << 85 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 109 // create G4Decayproducts << 86 //create G4Decayproducts 110 auto products = new G4DecayProducts(*parentp << 87 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 111 delete parentparticle; 88 delete parentparticle; 112 89 113 // calculate electron energy << 90 // calcurate electron energy 114 91 115 G4double michel_rho = 0.75; // Standard Mod << 92 G4double michel_rho = 0.75; //Standard Model Michel rho 116 G4double michel_delta = 0.75; // Standard M << 93 G4double michel_delta = 0.75; //Standard Model Michel delta 117 G4double michel_xsi = 1.00; // Standard Mod << 94 G4double michel_xsi = 1.00; //Standard Model Michel xsi 118 G4double michel_eta = 0.00; // Standard Mod << 95 G4double michel_eta = 0.00; //Standard Model eta 119 96 120 G4double rndm, x, ctheta; 97 G4double rndm, x, ctheta; 121 98 122 G4double FG; << 99 G4double FG; 123 G4double FG_max = 2.00; 100 G4double FG_max = 2.00; 124 101 125 G4double W_mue = (EMMU * EMMU + EMASS * EMAS << 102 G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU); 126 G4double x0 = EMASS / W_mue; << 103 G4double x0 = EMASS/W_mue; 127 104 128 G4double x0_squared = x0 * x0; << 105 G4double x0_squared = x0*x0; 129 106 130 // ***************************************** 107 // *************************************************** 131 // x0 <= x <= 1. and -1 <= y <= 1 108 // x0 <= x <= 1. and -1 <= y <= 1 132 // 109 // 133 // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g 110 // F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y 134 // ***************************************** 111 // *************************************************** 135 112 136 // ***** sampling F(x,y) directly (brute for 113 // ***** sampling F(x,y) directly (brute force) ***** 137 114 138 const std::size_t MAX_LOOP = 10000; << 115 do{ 139 for (std::size_t loop_count = 0; loop_count << 116 140 // Sample the positron energy by sampling 117 // Sample the positron energy by sampling from F 141 118 142 rndm = G4UniformRand(); 119 rndm = G4UniformRand(); 143 120 144 x = x0 + rndm * (1. - x0); << 121 x = x0 + rndm*(1.-x0); 145 122 146 G4double x_squared = x * x; << 123 G4double x_squared = x*x; 147 124 148 G4double F_IS, F_AS, G_IS, G_AS; 125 G4double F_IS, F_AS, G_IS, G_AS; 149 126 150 F_IS = 1. / 6. * (-2. * x_squared + 3. * x << 127 F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared); 151 F_AS = 1. / 6. * std::sqrt(x_squared - x0_ << 128 F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared)); 152 129 153 G_IS = 2. / 9. * (michel_rho - 0.75) * (4. << 130 G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared); 154 G_IS = G_IS + michel_eta * (1. - x) * x0; << 131 G_IS = G_IS + michel_eta*(1.-x)*x0; 155 132 156 G_AS = 3. * (michel_xsi - 1.) * (1. - x); << 133 G_AS = 3.*(michel_xsi-1.)*(1.-x); 157 G_AS = << 134 G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared)); 158 G_AS + 2. * (michel_xsi * michel_delta - << 135 G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS; 159 G_AS = 1. / 9. * std::sqrt(x_squared - x0_ << 160 136 161 F_IS = F_IS + G_IS; 137 F_IS = F_IS + G_IS; 162 F_AS = F_AS + G_AS; 138 F_AS = F_AS + G_AS; 163 139 164 // *** Radiative Corrections *** 140 // *** Radiative Corrections *** 165 141 166 const G4double omega = std::log(EMMU / EMA << 142 G4double R_IS = F_c(x,x0); 167 G4double R_IS = F_c(x, x0, omega); << 168 143 169 G4double F = 6. * F_IS + R_IS / std::sqrt( << 144 G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared); 170 145 171 // *** Radiative Corrections *** 146 // *** Radiative Corrections *** 172 147 173 G4double R_AS = F_theta(x, x0, omega); << 148 G4double R_AS = F_theta(x,x0); 174 149 175 rndm = G4UniformRand(); 150 rndm = G4UniformRand(); 176 151 177 ctheta = 2. * rndm - 1.; << 152 ctheta = 2.*rndm-1.; 178 153 179 G4double G = 6. * F_AS - R_AS / std::sqrt( << 154 G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared); 180 155 181 FG = std::sqrt(x_squared - x0_squared) * F << 156 FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta); 182 157 183 if (FG > FG_max) { << 158 if(FG>FG_max){ 184 G4Exception("G4MuonDecayChannelWithSpin: << 159 G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl; 185 "Problem in Muon Decay: FG > << 186 FG_max = FG; 160 FG_max = FG; 187 } 161 } 188 162 189 rndm = G4UniformRand(); 163 rndm = G4UniformRand(); 190 164 191 if (FG >= rndm * FG_max) break; << 165 }while(FG<rndm*FG_max); 192 } << 193 166 194 G4double energy = x * W_mue; 167 G4double energy = x * W_mue; 195 168 196 rndm = G4UniformRand(); 169 rndm = G4UniformRand(); 197 170 198 G4double phi = twopi * rndm; 171 G4double phi = twopi * rndm; 199 172 200 if (energy < EMASS) energy = EMASS; << 173 if(energy < EMASS) energy = EMASS; 201 174 202 // Calculate daughter momentum << 175 // calculate daughter momentum 203 G4double daughtermomentum[3]; 176 G4double daughtermomentum[3]; 204 177 205 daughtermomentum[0] = std::sqrt(energy * ene << 178 daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS); 206 179 207 G4double stheta = std::sqrt(1. - ctheta * ct << 180 G4double stheta = std::sqrt(1.-ctheta*ctheta); 208 G4double cphi = std::cos(phi); 181 G4double cphi = std::cos(phi); 209 G4double sphi = std::sin(phi); 182 G4double sphi = std::sin(phi); 210 183 211 // Coordinates of the decay positron with re << 184 //Coordinates of the decay positron with respect to the muon spin 212 G4double px = stheta * cphi; << 185 213 G4double py = stheta * sphi; << 186 G4double px = stheta*cphi; >> 187 G4double py = stheta*sphi; 214 G4double pz = ctheta; 188 G4double pz = ctheta; 215 189 216 G4ThreeVector direction0(px, py, pz); << 190 G4ThreeVector direction0(px,py,pz); 217 191 218 direction0.rotateUz(parent_polarization); 192 direction0.rotateUz(parent_polarization); 219 193 220 auto daughterparticle0 = << 194 G4DynamicParticle * daughterparticle0 221 new G4DynamicParticle(G4MT_daughters[0], d << 195 = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0); 222 196 223 products->PushProducts(daughterparticle0); 197 products->PushProducts(daughterparticle0); 224 198 >> 199 225 // daughter 1 ,2 (neutrinos) 200 // daughter 1 ,2 (neutrinos) 226 // create neutrinos in the C.M frame of two 201 // create neutrinos in the C.M frame of two neutrinos 227 G4double energy2 = parentmass - energy; << 202 G4double energy2 = parentmass*(1.0 - x/2.0); 228 G4double vmass = std::sqrt((energy2 - daught << 203 G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0])); 229 G4double beta = -1.0 * daughtermomentum[0] / << 204 G4double beta = -1.0*daughtermomentum[0]/energy2; 230 G4double costhetan = 2. * G4UniformRand() - << 205 G4double costhetan = 2.*G4UniformRand()-1.0; 231 G4double sinthetan = std::sqrt((1.0 - costhe << 206 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); 232 G4double phin = twopi * G4UniformRand() * ra << 207 G4double phin = twopi*G4UniformRand()*rad; 233 G4double sinphin = std::sin(phin); 208 G4double sinphin = std::sin(phin); 234 G4double cosphin = std::cos(phin); 209 G4double cosphin = std::cos(phin); 235 210 236 G4ThreeVector direction1(sinthetan * cosphin << 211 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan); 237 auto daughterparticle1 = new G4DynamicPartic << 212 G4DynamicParticle * daughterparticle1 238 auto daughterparticle2 = << 213 = new G4DynamicParticle( daughters[1], direction1*(vmass/2.)); 239 new G4DynamicParticle(G4MT_daughters[2], d << 214 G4DynamicParticle * daughterparticle2 >> 215 = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.)); 240 216 241 // boost to the muon rest frame 217 // boost to the muon rest frame 242 G4LorentzVector p4; 218 G4LorentzVector p4; 243 p4 = daughterparticle1->Get4Momentum(); 219 p4 = daughterparticle1->Get4Momentum(); 244 p4.boost(direction0.x() * beta, direction0.y << 220 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); 245 daughterparticle1->Set4Momentum(p4); 221 daughterparticle1->Set4Momentum(p4); 246 p4 = daughterparticle2->Get4Momentum(); 222 p4 = daughterparticle2->Get4Momentum(); 247 p4.boost(direction0.x() * beta, direction0.y << 223 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); 248 daughterparticle2->Set4Momentum(p4); 224 daughterparticle2->Set4Momentum(p4); 249 products->PushProducts(daughterparticle1); 225 products->PushProducts(daughterparticle1); 250 products->PushProducts(daughterparticle2); 226 products->PushProducts(daughterparticle2); 251 daughtermomentum[1] = daughterparticle1->Get 227 daughtermomentum[1] = daughterparticle1->GetTotalMomentum(); 252 daughtermomentum[2] = daughterparticle2->Get 228 daughtermomentum[2] = daughterparticle2->GetTotalMomentum(); 253 229 254 // output message 230 // output message 255 #ifdef G4VERBOSE 231 #ifdef G4VERBOSE 256 if (GetVerboseLevel() > 1) { << 232 if (GetVerboseLevel()>1) { 257 G4cout << "G4MuonDecayChannelWithSpin::Dec 233 G4cout << "G4MuonDecayChannelWithSpin::DecayIt "; 258 G4cout << " create decay products in rest << 234 G4cout << " create decay products in rest frame " <<G4endl; 259 G4double TT = daughterparticle0->GetTotalE << 235 products->DumpInfo(); 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 } 236 } 269 #endif 237 #endif 270 << 271 return products; 238 return products; 272 } 239 } 273 240 274 G4double G4MuonDecayChannelWithSpin::R_c(G4dou << 241 G4double G4MuonDecayChannelWithSpin::R_c(G4double x){ 275 { << 276 auto n_max = (G4int)(100. * x); << 277 242 278 if (n_max < 10) n_max = 10; << 243 G4int n_max = (int)(100.*x); >> 244 >> 245 if(n_max<10)n_max=10; 279 246 280 G4double L2 = 0.0; 247 G4double L2 = 0.0; 281 248 282 for (G4int n = 1; n <= n_max; ++n) { << 249 for(G4int n=1; n<=n_max; n++){ 283 L2 += std::pow(x, n) / (n * n); << 250 L2 += std::pow(x,n)/(n*n); 284 } 251 } 285 252 >> 253 G4double omega = std::log(EMMU/EMASS); >> 254 286 G4double r_c; 255 G4double r_c; 287 256 288 r_c = 2. * L2 - (pi * pi / 3.) - 2.; << 257 r_c = 2.*L2-(pi*pi/3.)-2.; 289 r_c = r_c + omega * (1.5 + 2. * std::log((1. << 258 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) << 259 r_c = r_c - std::log(x)*(2.*std::log(x)-1.); 291 r_c = r_c + (3. * std::log(x) - 1. - 1. / x) << 260 r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x); 292 261 293 return r_c; 262 return r_c; 294 } 263 } 295 264