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 // G4TauLeptonicDecayChannel class implementat 26 // G4TauLeptonicDecayChannel class implementation 27 // 27 // 28 // Author: H.Kurashige, 30 May 1997 << 28 // Author: H.Kurashige, 30 May 1997 29 // ------------------------------------------- 29 // -------------------------------------------------------------------- 30 30 31 #include "G4TauLeptonicDecayChannel.hh" << 32 << 33 #include "G4DecayProducts.hh" << 34 #include "G4LorentzRotation.hh" << 35 #include "G4LorentzVector.hh" << 36 #include "G4ParticleDefinition.hh" 31 #include "G4ParticleDefinition.hh" 37 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh" 38 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh" >> 34 #include "G4DecayProducts.hh" 39 #include "G4VDecayChannel.hh" 35 #include "G4VDecayChannel.hh" >> 36 #include "G4TauLeptonicDecayChannel.hh" 40 #include "Randomize.hh" 37 #include "Randomize.hh" >> 38 #include "G4LorentzVector.hh" >> 39 #include "G4LorentzRotation.hh" >> 40 41 41 42 G4TauLeptonicDecayChannel::G4TauLeptonicDecayC << 42 G4TauLeptonicDecayChannel::G4TauLeptonicDecayChannel() 43 << 43 : G4VDecayChannel() >> 44 { >> 45 } >> 46 >> 47 // -------------------------------------------------------------------- >> 48 G4TauLeptonicDecayChannel:: >> 49 G4TauLeptonicDecayChannel(const G4String& theParentName, >> 50 G4double theBR, >> 51 const G4String& theLeptonName) 44 : G4VDecayChannel("Tau Leptonic Decay", 1) 52 : G4VDecayChannel("Tau Leptonic Decay", 1) 45 { 53 { 46 // set names for daughter particles 54 // set names for daughter particles 47 if (theParentName == "tau+") { << 55 if (theParentName == "tau+") >> 56 { 48 SetBR(theBR); 57 SetBR(theBR); 49 SetParent("tau+"); 58 SetParent("tau+"); 50 SetNumberOfDaughters(3); 59 SetNumberOfDaughters(3); 51 if ((theLeptonName == "e-" || theLeptonNam << 60 if ((theLeptonName=="e-"||theLeptonName=="e+")) >> 61 { 52 SetDaughter(0, "e+"); 62 SetDaughter(0, "e+"); 53 SetDaughter(1, "nu_e"); 63 SetDaughter(1, "nu_e"); 54 SetDaughter(2, "anti_nu_tau"); 64 SetDaughter(2, "anti_nu_tau"); 55 } 65 } 56 else { << 66 else >> 67 { 57 SetDaughter(0, "mu+"); 68 SetDaughter(0, "mu+"); 58 SetDaughter(1, "nu_mu"); 69 SetDaughter(1, "nu_mu"); 59 SetDaughter(2, "anti_nu_tau"); 70 SetDaughter(2, "anti_nu_tau"); 60 } << 71 } 61 } 72 } 62 else if (theParentName == "tau-") { << 73 else if (theParentName == "tau-") >> 74 { 63 SetBR(theBR); 75 SetBR(theBR); 64 SetParent("tau-"); 76 SetParent("tau-"); 65 SetNumberOfDaughters(3); 77 SetNumberOfDaughters(3); 66 if ((theLeptonName == "e-" || theLeptonNam << 78 if ((theLeptonName=="e-"||theLeptonName=="e+")) >> 79 { 67 SetDaughter(0, "e-"); 80 SetDaughter(0, "e-"); 68 SetDaughter(1, "anti_nu_e"); 81 SetDaughter(1, "anti_nu_e"); 69 SetDaughter(2, "nu_tau"); 82 SetDaughter(2, "nu_tau"); 70 } 83 } 71 else { << 84 else >> 85 { 72 SetDaughter(0, "mu-"); 86 SetDaughter(0, "mu-"); 73 SetDaughter(1, "anti_nu_mu"); 87 SetDaughter(1, "anti_nu_mu"); 74 SetDaughter(2, "nu_tau"); 88 SetDaughter(2, "nu_tau"); 75 } << 89 } 76 } 90 } 77 else { << 91 else >> 92 { 78 #ifdef G4VERBOSE 93 #ifdef G4VERBOSE 79 if (GetVerboseLevel() > 0) { << 94 if (GetVerboseLevel()>0) >> 95 { 80 G4cout << "G4TauLeptonicDecayChannel:: c 96 G4cout << "G4TauLeptonicDecayChannel:: constructor :"; 81 G4cout << " parent particle is not tau b 97 G4cout << " parent particle is not tau but "; 82 G4cout << theParentName << G4endl; 98 G4cout << theParentName << G4endl; 83 } 99 } 84 #endif 100 #endif 85 } 101 } 86 } 102 } 87 103 >> 104 // -------------------------------------------------------------------- >> 105 G4TauLeptonicDecayChannel::~G4TauLeptonicDecayChannel() >> 106 { >> 107 } >> 108 >> 109 // -------------------------------------------------------------------- >> 110 G4TauLeptonicDecayChannel:: >> 111 G4TauLeptonicDecayChannel(const G4TauLeptonicDecayChannel& right) >> 112 : G4VDecayChannel(right) >> 113 { >> 114 } >> 115 >> 116 // -------------------------------------------------------------------- 88 G4TauLeptonicDecayChannel& 117 G4TauLeptonicDecayChannel& 89 G4TauLeptonicDecayChannel::operator=(const G4T 118 G4TauLeptonicDecayChannel::operator=(const G4TauLeptonicDecayChannel& right) 90 { 119 { 91 if (this != &right) { << 120 if (this != &right) >> 121 { 92 kinematics_name = right.kinematics_name; 122 kinematics_name = right.kinematics_name; 93 verboseLevel = right.verboseLevel; 123 verboseLevel = right.verboseLevel; 94 rbranch = right.rbranch; 124 rbranch = right.rbranch; 95 125 96 // copy parent name 126 // copy parent name 97 parent_name = new G4String(*right.parent_n 127 parent_name = new G4String(*right.parent_name); 98 128 99 // clear daughters_name array 129 // clear daughters_name array 100 ClearDaughtersName(); 130 ClearDaughtersName(); 101 131 102 // recreate array 132 // recreate array 103 numberOfDaughters = right.numberOfDaughter 133 numberOfDaughters = right.numberOfDaughters; 104 if (numberOfDaughters > 0) { << 134 if ( numberOfDaughters >0 ) 105 if (daughters_name != nullptr) ClearDaug << 135 { >> 136 if (daughters_name !=0) ClearDaughtersName(); 106 daughters_name = new G4String*[numberOfD 137 daughters_name = new G4String*[numberOfDaughters]; 107 // copy daughters name 138 // copy daughters name 108 for (G4int index = 0; index < numberOfDa << 139 for (G4int index=0; index<numberOfDaughters; ++index) >> 140 { 109 daughters_name[index] = new G4String(* 141 daughters_name[index] = new G4String(*right.daughters_name[index]); 110 } 142 } 111 } 143 } 112 } 144 } 113 return *this; 145 return *this; 114 } 146 } 115 147 116 G4DecayProducts* G4TauLeptonicDecayChannel::De << 148 // -------------------------------------------------------------------- >> 149 G4DecayProducts* G4TauLeptonicDecayChannel::DecayIt(G4double) 117 { 150 { 118 // this version neglects muon polarization << 151 // this version neglects muon polarization 119 // assumes the pure V-A couplin 152 // assumes the pure V-A coupling 120 // gives incorrect energy spect 153 // gives incorrect energy spectrum for neutrinos 121 154 122 #ifdef G4VERBOSE 155 #ifdef G4VERBOSE 123 if (GetVerboseLevel() > 1) G4cout << "G4TauL << 156 if (GetVerboseLevel()>1) G4cout << "G4TauLeptonicDecayChannel::DecayIt()"; 124 #endif 157 #endif 125 158 126 CheckAndFillParent(); 159 CheckAndFillParent(); 127 CheckAndFillDaughters(); 160 CheckAndFillDaughters(); 128 << 161 129 // parent mass 162 // parent mass 130 G4double parentmass = G4MT_parent->GetPDGMas 163 G4double parentmass = G4MT_parent->GetPDGMass(); 131 164 132 // daughters'mass 165 // daughters'mass 133 const G4int N_DAUGHTER = 3; << 166 const G4int N_DAUGHTER=3; 134 G4double daughtermass[N_DAUGHTER]; << 167 G4double daughtermass[N_DAUGHTER]; 135 for (G4int index = 0; index < N_DAUGHTER; ++ << 168 for (G4int index=0; index<N_DAUGHTER; ++index) >> 169 { 136 daughtermass[index] = G4MT_daughters[index 170 daughtermass[index] = G4MT_daughters[index]->GetPDGMass(); 137 } 171 } 138 172 139 // create parent G4DynamicParticle at rest 173 // create parent G4DynamicParticle at rest 140 G4ThreeVector dummy; 174 G4ThreeVector dummy; 141 auto parentparticle = new G4DynamicParticle( << 175 G4DynamicParticle* parentparticle >> 176 = new G4DynamicParticle(G4MT_parent, dummy, 0.0); 142 // create G4Decayproducts 177 // create G4Decayproducts 143 auto products = new G4DecayProducts(*parentp << 178 G4DecayProducts* products = new G4DecayProducts(*parentparticle); 144 delete parentparticle; 179 delete parentparticle; 145 180 146 // calculate daughter momentum 181 // calculate daughter momentum 147 G4double daughtermomentum[N_DAUGHTER]; 182 G4double daughtermomentum[N_DAUGHTER]; 148 183 149 // calculate lepton momentum 184 // calculate lepton momentum 150 G4double pmax = (parentmass * parentmass - d << 185 G4double pmax = (parentmass*parentmass >> 186 - daughtermass[0]*daughtermass[0])/2./parentmass; 151 G4double p, e; 187 G4double p, e; 152 G4double r; 188 G4double r; 153 const std::size_t MAX_LOOP = 10000; << 189 const std::size_t MAX_LOOP=10000; 154 for (std::size_t loop_counter = 0; loop_coun << 190 for (std::size_t loop_counter=0; loop_counter<MAX_LOOP; ++loop_counter) >> 191 { 155 // determine momentum/energy 192 // determine momentum/energy 156 r = G4UniformRand(); 193 r = G4UniformRand(); 157 p = pmax * G4UniformRand(); << 194 p = pmax*G4UniformRand(); 158 e = std::sqrt(p * p + daughtermass[0] * da << 195 e = std::sqrt(p*p + daughtermass[0]*daughtermass[0]); 159 if (r < spectrum(p, e, parentmass, daughte << 196 if (r < spectrum(p,e,parentmass,daughtermass[0]) ) break; 160 } << 197 } 161 198 162 // create daughter G4DynamicParticle << 199 // create daughter G4DynamicParticle 163 // daughter 0 (lepton) 200 // daughter 0 (lepton) 164 daughtermomentum[0] = p; 201 daughtermomentum[0] = p; 165 G4double costheta, sintheta, phi, sinphi, co << 202 G4double costheta, sintheta, phi, sinphi, cosphi; 166 costheta = 2. * G4UniformRand() - 1.0; << 203 costheta = 2.*G4UniformRand()-1.0; 167 sintheta = std::sqrt((1.0 - costheta) * (1.0 << 204 sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); 168 phi = twopi * G4UniformRand() * rad; << 205 phi = twopi*G4UniformRand()*rad; 169 sinphi = std::sin(phi); 206 sinphi = std::sin(phi); 170 cosphi = std::cos(phi); 207 cosphi = std::cos(phi); 171 G4ThreeVector direction0(sintheta * cosphi, << 208 G4ThreeVector direction0(sintheta*cosphi,sintheta*sinphi,costheta); 172 auto daughterparticle = << 209 G4DynamicParticle * daughterparticle 173 new G4DynamicParticle(G4MT_daughters[0], d << 210 = new G4DynamicParticle(G4MT_daughters[0], direction0*daughtermomentum[0]); 174 products->PushProducts(daughterparticle); 211 products->PushProducts(daughterparticle); 175 212 176 // daughter 1 ,2 (nutrinos) 213 // daughter 1 ,2 (nutrinos) 177 // create neutrinos in the C.M frame of two 214 // create neutrinos in the C.M frame of two neutrinos 178 G4double energy2 = parentmass - e; << 215 G4double energy2 = parentmass-e; 179 G4double vmass = std::sqrt((energy2 - daught << 216 G4double vmass = std::sqrt((energy2-daughtermomentum[0]) 180 G4double beta = -1.0 * daughtermomentum[0] / << 217 *(energy2+daughtermomentum[0])); 181 G4double costhetan = 2. * G4UniformRand() - << 218 G4double beta = -1.0*daughtermomentum[0]/energy2; 182 G4double sinthetan = std::sqrt((1.0 - costhe << 219 G4double costhetan = 2.*G4UniformRand()-1.0; 183 G4double phin = twopi * G4UniformRand() * ra << 220 G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan)); >> 221 G4double phin = twopi*G4UniformRand()*rad; 184 G4double sinphin = std::sin(phin); 222 G4double sinphin = std::sin(phin); 185 G4double cosphin = std::cos(phin); 223 G4double cosphin = std::cos(phin); 186 224 187 G4ThreeVector direction1(sinthetan * cosphin << 225 G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan); 188 auto daughterparticle1 = new G4DynamicPartic << 226 G4DynamicParticle * daughterparticle1 189 auto daughterparticle2 = << 227 = new G4DynamicParticle( G4MT_daughters[1], direction1*(vmass/2.)); 190 new G4DynamicParticle(G4MT_daughters[2], d << 228 G4DynamicParticle * daughterparticle2 >> 229 = new G4DynamicParticle( G4MT_daughters[2], direction1*(-1.0*vmass/2.)); 191 230 192 // boost to the muon rest frame 231 // boost to the muon rest frame 193 G4LorentzVector p4; 232 G4LorentzVector p4; 194 p4 = daughterparticle1->Get4Momentum(); 233 p4 = daughterparticle1->Get4Momentum(); 195 p4.boost(direction0.x() * beta, direction0.y << 234 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); 196 daughterparticle1->Set4Momentum(p4); 235 daughterparticle1->Set4Momentum(p4); 197 p4 = daughterparticle2->Get4Momentum(); 236 p4 = daughterparticle2->Get4Momentum(); 198 p4.boost(direction0.x() * beta, direction0.y << 237 p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta); 199 daughterparticle2->Set4Momentum(p4); 238 daughterparticle2->Set4Momentum(p4); 200 products->PushProducts(daughterparticle1); 239 products->PushProducts(daughterparticle1); 201 products->PushProducts(daughterparticle2); 240 products->PushProducts(daughterparticle2); 202 daughtermomentum[1] = daughterparticle1->Get 241 daughtermomentum[1] = daughterparticle1->GetTotalMomentum(); 203 daughtermomentum[2] = daughterparticle2->Get 242 daughtermomentum[2] = daughterparticle2->GetTotalMomentum(); 204 243 205 // output message 244 // output message 206 #ifdef G4VERBOSE 245 #ifdef G4VERBOSE 207 if (GetVerboseLevel() > 1) { << 246 if (GetVerboseLevel()>1) >> 247 { 208 G4cout << "G4TauLeptonicDecayChannel::Deca 248 G4cout << "G4TauLeptonicDecayChannel::DecayIt "; 209 G4cout << " create decay products in rest 249 G4cout << " create decay products in rest frame " << G4endl; 210 products->DumpInfo(); 250 products->DumpInfo(); 211 } 251 } 212 #endif 252 #endif 213 return products; 253 return products; 214 } 254 } 215 255 216 G4double G4TauLeptonicDecayChannel::spectrum(G << 256 // -------------------------------------------------------------------- >> 257 G4double G4TauLeptonicDecayChannel::spectrum(G4double p, >> 258 G4double e, >> 259 G4double mtau, >> 260 G4double ml) 217 { 261 { 218 G4double f1; 262 G4double f1; 219 f1 = 3.0 * e * (mtau * mtau + ml * ml) - 4.0 << 263 f1 = 3.0*e*(mtau*mtau+ml*ml)-4.0*mtau*e*e-2.0*mtau*ml*ml; 220 return p * (f1) / (mtau * mtau * mtau * mtau << 264 return p*(f1)/(mtau*mtau*mtau*mtau)/(0.6); 221 } 265 } 222 266