Geant4 Cross Reference |
>> 1 // This code implementation is the intellectual property of >> 2 // the GEANT4 collaboration. 1 // 3 // 2 // ******************************************* << 4 // By copying, distributing or modifying the Program (or any work 3 // * License and Disclaimer << 5 // based on the Program) you indicate your acceptance of this statement, 4 // * << 6 // and all its terms. 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 // 7 // 26 // G4MuonDecayChannel class implementation << 8 // $Id: G4MuonDecayChannel.cc,v 1.4.8.1 1999/12/07 20:49:56 gunter Exp $ >> 9 // GEANT4 tag $Name: geant4-01-00 $ 27 // 10 // 28 // Author: H.Kurashige, 30 May 1997 << 11 // 29 // Contributions: << 12 // ------------------------------------------------------------ 30 // - 2005 - M.Melissas, J.Brunner CPPM/IN2P3 << 13 // GEANT 4 class header file 31 // Added V-A fluxes for neutrinos using a ne << 14 // 32 // ------------------------------------------- << 15 // For information related to this code contact: 33 << 16 // CERN, CN Division, ASD group 34 #include "G4MuonDecayChannel.hh" << 17 // History: first implementation, based on object model of >> 18 // 30 May 1997 H.Kurashige >> 19 // 10 June 1997 H.Kurashige >> 20 // ------------------------------------------------------------ 35 21 36 #include "G4DecayProducts.hh" << 37 #include "G4LorentzRotation.hh" << 38 #include "G4LorentzVector.hh" << 39 #include "G4ParticleDefinition.hh" 22 #include "G4ParticleDefinition.hh" 40 #include "G4PhysicalConstants.hh" << 23 #include "G4DecayProducts.hh" 41 #include "G4RotationMatrix.hh" << 42 #include "G4SystemOfUnits.hh" << 43 #include "G4VDecayChannel.hh" 24 #include "G4VDecayChannel.hh" >> 25 #include "G4MuonDecayChannel.hh" 44 #include "Randomize.hh" 26 #include "Randomize.hh" >> 27 #include "G4LorentzVector.hh" >> 28 #include "G4LorentzRotation.hh" >> 29 45 30 46 G4MuonDecayChannel::G4MuonDecayChannel(const G << 31 G4MuonDecayChannel::G4MuonDecayChannel(const G4String& theParentName, 47 : G4VDecayChannel("Muon Decay", 1) << 32 G4double theBR) >> 33 :G4VDecayChannel("Muon Decay",1) 48 { 34 { >> 35 //#ifdef G4VERBOSE >> 36 //if (GetVerboseLevel()>1) { >> 37 // G4cout << "G4MuonDecayChannel:: constructor "; >> 38 // G4cout << "addr[" << this << "]" << endl; >> 39 //} >> 40 //#endif >> 41 49 // set names for daughter particles 42 // set names for daughter particles 50 if (theParentName == "mu+") { 43 if (theParentName == "mu+") { 51 SetBR(theBR); 44 SetBR(theBR); 52 SetParent("mu+"); 45 SetParent("mu+"); 53 SetNumberOfDaughters(3); 46 SetNumberOfDaughters(3); 54 SetDaughter(0, "e+"); 47 SetDaughter(0, "e+"); 55 SetDaughter(1, "nu_e"); 48 SetDaughter(1, "nu_e"); 56 SetDaughter(2, "anti_nu_mu"); 49 SetDaughter(2, "anti_nu_mu"); 57 } << 50 } else if (theParentName == "mu-") { 58 else if (theParentName == "mu-") { << 59 SetBR(theBR); 51 SetBR(theBR); 60 SetParent("mu-"); 52 SetParent("mu-"); 61 SetNumberOfDaughters(3); 53 SetNumberOfDaughters(3); 62 SetDaughter(0, "e-"); 54 SetDaughter(0, "e-"); 63 SetDaughter(1, "anti_nu_e"); 55 SetDaughter(1, "anti_nu_e"); 64 SetDaughter(2, "nu_mu"); 56 SetDaughter(2, "nu_mu"); 65 } << 57 } else { 66 else { << 58 //#ifdef G4VERBOSE 67 #ifdef G4VERBOSE << 59 // if (GetVerboseLevel()>0) { 68 if (GetVerboseLevel() > 0) { << 60 // G4cout << "G4MuonDecayChannel:: constructor :"; 69 G4cout << "G4MuonDecayChannel:: construc << 61 // G4cout << " parent particle is not muon but "; 70 G4cout << " parent particle is not muon << 62 // G4cout << theParentName << endl; 71 G4cout << theParentName << G4endl; << 63 // } 72 } << 64 //#endif 73 #endif << 74 } 65 } 75 } 66 } 76 67 77 G4MuonDecayChannel& G4MuonDecayChannel::operat << 68 G4MuonDecayChannel::~G4MuonDecayChannel() 78 { 69 { 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 } 70 } 103 71 104 G4DecayProducts* G4MuonDecayChannel::DecayIt(G << 72 G4DecayProducts *G4MuonDecayChannel::DecayIt(G4double) 105 { 73 { 106 // this version neglects muon polarization,a << 74 // this version neglects muon polarization 107 // assumes the pure V-A couplin 75 // assumes the pure V-A coupling 108 // the Neutrinos are correctly << 76 // gives incorrect energy spectrum for neutrinos 109 << 110 #ifdef G4VERBOSE 77 #ifdef G4VERBOSE 111 if (GetVerboseLevel() > 1) G4cout << "G4Muon << 78 if (GetVerboseLevel()>1) G4cout << "G4MuonDecayChannel::DecayIt "; 112 #endif 79 #endif 113 80 114 CheckAndFillParent(); << 81 if (parent == 0) FillParent(); 115 CheckAndFillDaughters(); << 82 if (daughters == 0) FillDaughters(); 116 << 83 117 // parent mass 84 // parent mass 118 G4double parentmass = G4MT_parent->GetPDGMas << 85 G4double parentmass = parent->GetPDGMass(); 119 const G4int N_DAUGHTER = 3; << 120 86 121 // daughters'mass << 87 //daughters'mass 122 G4double daughtermass[N_DAUGHTER]; << 88 G4double daughtermass[3]; 123 // G4double sumofdaughtermass = 0.0; << 89 G4double sumofdaughtermass = 0.0; 124 for (G4int index = 0; index < N_DAUGHTER; ++ << 90 for (G4int index=0; index<3; index++){ 125 daughtermass[index] = G4MT_daughters[index << 91 daughtermass[index] = daughters[index]->GetPDGMass(); 126 // sumofdaughtermass += daughtermass[index << 92 sumofdaughtermass += daughtermass[index]; 127 } 93 } 128 94 129 // create parent G4DynamicParticle at rest << 95 //create parent G4DynamicParticle at rest 130 G4ThreeVector dummy; 96 G4ThreeVector dummy; 131 auto parentparticle = new G4DynamicParticle( << 97 G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0); 132 // create G4Decayproducts << 98 //create G4Decayproducts 133 auto products = new G4DecayProducts(*parentp << 99 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 134 delete parentparticle; 100 delete parentparticle; 135 101 136 // calculate daughter momentum 102 // calculate daughter momentum 137 G4double daughtermomentum[N_DAUGHTER]; << 103 G4double daughtermomentum[3]; 138 // calculate electron energy << 104 G4double momentumsum = 0.0; 139 G4double xmax = (1.0 + daughtermass[0] * dau << 105 G4double energy; >> 106 // calcurate electron energy >> 107 G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass); 140 G4double x; 108 G4double x; 141 << 109 G4double r; 142 G4double Ee, Ene; << 110 do { 143 << 111 do { 144 G4double gam; << 112 r = G4UniformRand(); 145 G4double EMax = parentmass / 2 - daughtermas << 113 x = xmax*G4UniformRand(); 146 << 114 } while (r < (3.0 - 2.0*x)*x*x); 147 const std::size_t MAX_LOOP = 1000; << 115 energy = x*parentmass/2.0 - daughtermass[0]; 148 // Generating Random Energy << 116 } while (energy <0.0); 149 for (std::size_t loop1 = 0; loop1 < MAX_LOOP << 117 //create daughter G4DynamicParticle 150 Ee = G4UniformRand(); << 118 // daughter 0 (electron) 151 for (std::size_t loop2 = 0; loop2 < MAX_LO << 119 daughtermomentum[0] = sqrt(energy*energy + 2.0*energy* daughtermass[0]); 152 x = xmax * G4UniformRand(); << 120 G4double costheta, sintheta, phi, sinphi, cosphi; 153 gam = G4UniformRand(); << 121 costheta = 2.*G4UniformRand()-1.0; 154 if (gam <= x * (1. - x)) break; << 122 sintheta = sqrt((1.0-costheta)*(1.0+costheta)); 155 x = xmax; << 123 phi = 2.0*M_PI*G4UniformRand()*rad; 156 } << 124 sinphi = sin(phi); 157 Ene = x; << 125 cosphi = cos(phi); 158 if (Ene >= (1. - Ee)) break; << 126 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta); 159 Ene = 1. - Ee; << 127 G4DynamicParticle * daughterparticle 160 } << 128 = new G4DynamicParticle( daughters[0], direction0*daughtermomentum[0]); 161 G4double Enm = (2. - Ee - Ene); << 162 << 163 // initialisation of rotation parameters << 164 << 165 G4double costheta, sintheta, rphi, rtheta, r << 166 costheta = 1. - 2. / Ee - 2. / Ene + 2. / En << 167 sintheta = std::sqrt(1. - costheta * costhet << 168 << 169 rphi = twopi * G4UniformRand() * rad; << 170 rtheta = (std::acos(2. * G4UniformRand() - 1 << 171 rpsi = twopi * G4UniformRand() * rad; << 172 << 173 G4RotationMatrix rot; << 174 rot.set(rphi, rtheta, rpsi); << 175 << 176 // electron 0 << 177 daughtermomentum[0] = std::sqrt(Ee * Ee * EM << 178 G4ThreeVector direction0(0.0, 0.0, 1.0); << 179 << 180 direction0 *= rot; << 181 << 182 auto daughterparticle = << 183 new G4DynamicParticle(G4MT_daughters[0], d << 184 << 185 products->PushProducts(daughterparticle); 129 products->PushProducts(daughterparticle); 186 130 187 // electronic neutrino 1 << 131 // daughter 1 ,2 (nutrinos) 188 << 132 // create neutrinos in the C.M frame of two neutrinos 189 daughtermomentum[1] = std::sqrt(Ene * Ene * << 133 G4double energy2 = parentmass*(1.0 - x/2.0); 190 G4ThreeVector direction1(sintheta, 0.0, cost << 134 G4double vmass = sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0])); 191 << 135 G4double beta = -1.0*daughtermomentum[0]/energy2; 192 direction1 *= rot; << 136 G4double costhetan = 2.*G4UniformRand()-1.0; 193 << 137 G4double sinthetan = sqrt((1.0-costhetan)*(1.0+costhetan)); 194 auto daughterparticle1 = << 138 G4double phin = 2.0*M_PI*G4UniformRand()*rad; 195 new G4DynamicParticle(G4MT_daughters[1], d << 139 G4double sinphin = sin(phin); >> 140 G4double cosphin = cos(phin); >> 141 >> 142 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan); >> 143 G4DynamicParticle * daughterparticle1 >> 144 = new G4DynamicParticle( daughters[1], direction1*(vmass/2.)); >> 145 G4DynamicParticle * daughterparticle2 >> 146 = new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.)); >> 147 >> 148 // boost to the muon rest frame >> 149 G4LorentzVector p4; >> 150 p4 = daughterparticle1->Get4Momentum(); >> 151 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); >> 152 daughterparticle1->Set4Momentum(p4); >> 153 p4 = daughterparticle2->Get4Momentum(); >> 154 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); >> 155 daughterparticle2->Set4Momentum(p4); 196 products->PushProducts(daughterparticle1); 156 products->PushProducts(daughterparticle1); 197 << 198 // muonnic neutrino 2 << 199 << 200 daughtermomentum[2] = std::sqrt(Enm * Enm * << 201 G4ThreeVector direction2(-Ene / Enm * sinthe << 202 << 203 direction2 *= rot; << 204 << 205 auto daughterparticle2 = << 206 new G4DynamicParticle(G4MT_daughters[2], d << 207 products->PushProducts(daughterparticle2); 157 products->PushProducts(daughterparticle2); >> 158 daughtermomentum[1] = daughterparticle1->GetTotalMomentum(); >> 159 daughtermomentum[2] = daughterparticle2->GetTotalMomentum(); 208 160 209 // output message << 161 // output message 210 #ifdef G4VERBOSE 162 #ifdef G4VERBOSE 211 if (GetVerboseLevel() > 1) { << 163 if (GetVerboseLevel()>1) { 212 G4cout << "G4MuonDecayChannel::DecayIt()"; << 164 G4cout << "G4MuonDecayChannel::DecayIt "; 213 G4cout << " create decay products in rest << 165 G4cout << " create decay products in rest frame " <<endl; 214 products->DumpInfo(); 166 products->DumpInfo(); 215 } 167 } 216 #endif 168 #endif 217 return products; 169 return products; 218 } 170 } >> 171 >> 172 >> 173 >> 174 >> 175 >> 176 219 177