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 // G4MuonDecayChannel class implementation << 27 // 26 // 28 // Author: H.Kurashige, 30 May 1997 << 27 // $Id: G4MuonDecayChannel.cc,v 1.17 2006/06/29 19:25:34 gunter Exp $ 29 // Contributions: << 28 // GEANT4 tag $Name: geant4-09-02 $ 30 // - 2005 - M.Melissas, J.Brunner CPPM/IN2P3 << 29 // 31 // Added V-A fluxes for neutrinos using a ne << 30 // 32 // ------------------------------------------- << 31 // ------------------------------------------------------------ 33 << 32 // GEANT 4 class header file 34 #include "G4MuonDecayChannel.hh" << 33 // >> 34 // History: first implementation, based on object model of >> 35 // 30 May 1997 H.Kurashige >> 36 // >> 37 // Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige >> 38 //2005 >> 39 // M. Melissas ( melissas AT cppm.in2p3.fr) >> 40 // J. Brunner ( brunner AT cppm.in2p3.fr) >> 41 // Adding V-A fluxes for neutrinos using a new algortithm : >> 42 // ------------------------------------------------------------ 35 43 36 #include "G4DecayProducts.hh" << 37 #include "G4LorentzRotation.hh" << 38 #include "G4LorentzVector.hh" << 39 #include "G4ParticleDefinition.hh" 44 #include "G4ParticleDefinition.hh" 40 #include "G4PhysicalConstants.hh" << 45 #include "G4DecayProducts.hh" 41 #include "G4RotationMatrix.hh" << 42 #include "G4SystemOfUnits.hh" << 43 #include "G4VDecayChannel.hh" 46 #include "G4VDecayChannel.hh" >> 47 #include "G4MuonDecayChannel.hh" 44 #include "Randomize.hh" 48 #include "Randomize.hh" >> 49 #include "G4LorentzVector.hh" >> 50 #include "G4LorentzRotation.hh" >> 51 #include "G4RotationMatrix.hh" >> 52 45 53 46 G4MuonDecayChannel::G4MuonDecayChannel(const G << 54 G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName, 47 : G4VDecayChannel("Muon Decay", 1) << 55 G4double theBR) >> 56 :G4VDecayChannel("Muon Decay",1) 48 { 57 { 49 // set names for daughter particles 58 // set names for daughter particles 50 if (theParentName == "mu+") { 59 if (theParentName == "mu+") { 51 SetBR(theBR); 60 SetBR(theBR); 52 SetParent("mu+"); 61 SetParent("mu+"); 53 SetNumberOfDaughters(3); 62 SetNumberOfDaughters(3); 54 SetDaughter(0, "e+"); 63 SetDaughter(0, "e+"); 55 SetDaughter(1, "nu_e"); 64 SetDaughter(1, "nu_e"); 56 SetDaughter(2, "anti_nu_mu"); 65 SetDaughter(2, "anti_nu_mu"); 57 } << 66 } else if (theParentName == "mu-") { 58 else if (theParentName == "mu-") { << 59 SetBR(theBR); 67 SetBR(theBR); 60 SetParent("mu-"); 68 SetParent("mu-"); 61 SetNumberOfDaughters(3); 69 SetNumberOfDaughters(3); 62 SetDaughter(0, "e-"); 70 SetDaughter(0, "e-"); 63 SetDaughter(1, "anti_nu_e"); 71 SetDaughter(1, "anti_nu_e"); 64 SetDaughter(2, "nu_mu"); 72 SetDaughter(2, "nu_mu"); 65 } << 73 } else { 66 else { << 67 #ifdef G4VERBOSE 74 #ifdef G4VERBOSE 68 if (GetVerboseLevel() > 0) { << 75 if (GetVerboseLevel()>0) { 69 G4cout << "G4MuonDecayChannel:: construc 76 G4cout << "G4MuonDecayChannel:: constructor :"; 70 G4cout << " parent particle is not muon 77 G4cout << " parent particle is not muon but "; 71 G4cout << theParentName << G4endl; 78 G4cout << theParentName << G4endl; 72 } 79 } 73 #endif 80 #endif 74 } 81 } 75 } 82 } 76 83 77 G4MuonDecayChannel& G4MuonDecayChannel::operat << 84 G4MuonDecayChannel::~G4MuonDecayChannel() 78 { 85 { 79 if (this != &right) { << 80 kinematics_name = right.kinematics_name; << 81 verboseLevel = right.verboseLevel; << 82 rbranch = right.rbranch; << 83 << 84 // copy parent name << 85 parent_name = new G4String(*right.parent_n << 86 << 87 // clear daughters_name array << 88 ClearDaughtersName(); << 89 << 90 // recreate array << 91 numberOfDaughters = right.numberOfDaughter << 92 if (numberOfDaughters > 0) { << 93 if (daughters_name != nullptr) ClearDaug << 94 daughters_name = new G4String*[numberOfD << 95 // copy daughters name << 96 for (G4int index = 0; index < numberOfDa << 97 daughters_name[index] = new G4String(* << 98 } << 99 } << 100 } << 101 return *this; << 102 } 86 } 103 87 104 G4DecayProducts* G4MuonDecayChannel::DecayIt(G << 88 G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double) 105 { 89 { 106 // this version neglects muon polarization,a << 90 // this version neglects muon polarization,and electron mass 107 // assumes the pure V-A couplin 91 // assumes the pure V-A coupling 108 // the Neutrinos are correctly << 92 // the Neutrinos are correctly V-A. 109 << 110 #ifdef G4VERBOSE 93 #ifdef G4VERBOSE 111 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 94 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt "; 112 #endif 95 #endif 113 96 114 CheckAndFillParent(); << 97 if (parent == 0) FillParent(); 115 CheckAndFillDaughters(); << 98 if (daughters == 0) FillDaughters(); 116 << 99 117 // parent mass 100 // parent mass 118 G4double parentmass = G4MT_parent->GetPDGMas << 101 G4double parentmass = parent->GetPDGMass(); 119 const G4int N_DAUGHTER = 3; << 120 102 121 // daughters'mass << 103 //daughters'mass 122 G4double daughtermass[N_DAUGHTER]; << 104 G4double daughtermass[3]; 123 // G4double sumofdaughtermass = 0.0; << 105 G4double sumofdaughtermass = 0.0; 124 for (G4int index = 0; index < N_DAUGHTER; ++ << 106 for (G4int index=0; index<3; index++){ 125 daughtermass[index] = G4MT_daughters[index << 107 daughtermass[index] = daughters[index]->GetPDGMass(); 126 // sumofdaughtermass += daughtermass[index << 108 sumofdaughtermass += daughtermass[index]; 127 } 109 } 128 110 129 // create parent G4DynamicParticle at rest << 111 //create parent G4DynamicParticle at rest 130 G4ThreeVector dummy; 112 G4ThreeVector dummy; 131 auto parentparticle = new G4DynamicParticle( << 113 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 132 // create G4Decayproducts << 114 //create G4Decayproducts 133 auto products = new G4DecayProducts(*parentp << 115 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 134 delete parentparticle; 116 delete parentparticle; 135 117 136 // calculate daughter momentum 118 // calculate daughter momentum 137 G4double daughtermomentum[N_DAUGHTER]; << 119 G4double daughtermomentum[3]; 138 // calculate electron energy << 120 // calcurate electron energy 139 G4double xmax = (1.0 + daughtermass[0] * dau << 121 G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass); 140 G4double x; 122 G4double x; 141 << 123 142 G4double Ee, Ene; << 124 G4double Ee,Ene; 143 << 125 144 G4double gam; 126 G4double gam; 145 G4double EMax = parentmass / 2 - daughtermas << 127 G4double EMax=parentmass/2-daughtermass[0]; 146 << 128 147 const std::size_t MAX_LOOP = 1000; << 129 148 // Generating Random Energy << 130 //Generating Random Energy 149 for (std::size_t loop1 = 0; loop1 < MAX_LOOP << 131 do { 150 Ee = G4UniformRand(); << 132 Ee=G4UniformRand(); 151 for (std::size_t loop2 = 0; loop2 < MAX_LO << 133 do{ 152 x = xmax * G4UniformRand(); << 134 x=xmax*G4UniformRand(); 153 gam = G4UniformRand(); << 135 gam=G4UniformRand(); 154 if (gam <= x * (1. - x)) break; << 136 }while (gam >x*(1.-x)); 155 x = xmax; << 137 Ene=x; 156 } << 138 } while ( Ene < (1.-Ee)); 157 Ene = x; << 139 G4double Enm=(2.-Ee-Ene); 158 if (Ene >= (1. - Ee)) break; << 140 159 Ene = 1. - Ee; << 141 160 } << 142 //initialisation of rotation parameters 161 G4double Enm = (2. - Ee - Ene); << 143 162 << 144 G4double costheta,sintheta,rphi,rtheta,rpsi; 163 // initialisation of rotation parameters << 145 costheta= 1.-2./Ee-2./Ene+2./Ene/Ee; 164 << 146 sintheta=std::sqrt(1.-costheta*costheta); 165 G4double costheta, sintheta, rphi, rtheta, r << 147 166 costheta = 1. - 2. / Ee - 2. / Ene + 2. / En << 148 167 sintheta = std::sqrt(1. - costheta * costhet << 149 rphi=twopi*G4UniformRand()*rad; 168 << 150 rtheta=(std::acos(2.*G4UniformRand()-1.)); 169 rphi = twopi * G4UniformRand() * rad; << 151 rpsi=twopi*G4UniformRand()*rad; 170 rtheta = (std::acos(2. * G4UniformRand() - 1 << 171 rpsi = twopi * G4UniformRand() * rad; << 172 152 173 G4RotationMatrix rot; 153 G4RotationMatrix rot; 174 rot.set(rphi, rtheta, rpsi); << 154 rot.set(rphi,rtheta,rpsi); 175 155 176 // electron 0 << 156 //electron 0 177 daughtermomentum[0] = std::sqrt(Ee * Ee * EM << 157 daughtermomentum[0]=std::sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]); 178 G4ThreeVector direction0(0.0, 0.0, 1.0); << 158 G4ThreeVector direction0(0.0,0.0,1.0); 179 159 180 direction0 *= rot; 160 direction0 *= rot; 181 161 182 auto daughterparticle = << 162 G4DynamicParticle * daughterparticle = new G4DynamicParticle ( daughters[0], direction0 * daughtermomentum[0]); 183 new G4DynamicParticle(G4MT_daughters[0], d << 184 163 185 products->PushProducts(daughterparticle); 164 products->PushProducts(daughterparticle); >> 165 >> 166 //electronic neutrino 1 186 167 187 // electronic neutrino 1 << 168 daughtermomentum[1]=std::sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]); 188 << 169 G4ThreeVector direction1(sintheta,0.0,costheta); 189 daughtermomentum[1] = std::sqrt(Ene * Ene * << 190 G4ThreeVector direction1(sintheta, 0.0, cost << 191 170 192 direction1 *= rot; 171 direction1 *= rot; 193 172 194 auto daughterparticle1 = << 173 G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( daughters[1], direction1 * daughtermomentum[1]); 195 new G4DynamicParticle(G4MT_daughters[1], d << 196 products->PushProducts(daughterparticle1); 174 products->PushProducts(daughterparticle1); 197 175 198 // muonnic neutrino 2 << 176 //muonnic neutrino 2 199 << 177 200 daughtermomentum[2] = std::sqrt(Enm * Enm * << 178 daughtermomentum[2]=std::sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]); 201 G4ThreeVector direction2(-Ene / Enm * sinthe << 179 G4ThreeVector direction2(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta); 202 180 203 direction2 *= rot; 181 direction2 *= rot; 204 182 205 auto daughterparticle2 = << 183 G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( daughters[2], 206 new G4DynamicParticle(G4MT_daughters[2], d << 184 direction2 * daughtermomentum[2]); 207 products->PushProducts(daughterparticle2); 185 products->PushProducts(daughterparticle2); 208 186 209 // output message << 187 >> 188 >> 189 >> 190 // output message 210 #ifdef G4VERBOSE 191 #ifdef G4VERBOSE 211 if (GetVerboseLevel() > 1) { << 192 if (GetVerboseLevel()>1) { 212 G4cout << "G4MuonDecayChannel::DecayIt()"; << 193 G4cout << "G4MuonDecayChannel::DecayIt "; 213 G4cout << " create decay products in rest << 194 G4cout << " create decay products in rest frame " <<G4endl; 214 products->DumpInfo(); 195 products->DumpInfo(); 215 } 196 } 216 #endif 197 #endif 217 return products; 198 return products; 218 } 199 } >> 200 >> 201 >> 202 >> 203 >> 204 >> 205 219 206