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 // G4PionRadiativeDecayChannel class implement << 26 // ------------------------------------------------------------ 27 // GEANT 4 class header file 27 // GEANT 4 class header file 28 // 28 // 29 // Author: P.Gumplinger, 30 July 2007 << 29 // History: 30 // Reference: M. Blecher, TRIUMF/PIENU Technot << 30 // 01 August 2007 P.Gumplinger 31 // "Inclusion of pi->enug in the Mo << 31 // Reference: TRIUMF PIENU Technote: 32 // ------------------------------------------- << 32 // M. Blecher - "Inclusion of pi->enug in MC " >> 33 // Rate is for gammas > 100keV >> 34 // >> 35 // ------------------------------------------------------------ >> 36 // >> 37 // >> 38 // 33 39 34 #include "G4PionRadiativeDecayChannel.hh" 40 #include "G4PionRadiativeDecayChannel.hh" 35 41 >> 42 #include "Randomize.hh" 36 #include "G4DecayProducts.hh" 43 #include "G4DecayProducts.hh" 37 #include "G4LorentzVector.hh" 44 #include "G4LorentzVector.hh" 38 #include "G4PhysicalConstants.hh" << 39 #include "G4SystemOfUnits.hh" << 40 #include "Randomize.hh" << 41 << 42 namespace << 43 { << 44 const G4double beta = 3.6612e-03; << 45 const G4double cib = 1.16141e-03; << 46 const G4double csdp = 3.45055e-02; << 47 const G4double csdm = 5.14122e-03; << 48 const G4double cif = 4.63543e-05; << 49 const G4double cig = 1.78928e-05; << 50 const G4double xl = 2. * 0.1 * MeV / 139.57 * << 51 const G4double yl = ((1. - xl) + std::sqrt((1 << 52 << 53 const G4double xu = 1. - (yl - std::sqrt(yl * << 54 const G4double yu = 1. + beta * beta; << 55 << 56 inline G4double D2W(const G4double x, const G4 << 57 { << 58 return cib * (1. - y) * (1. + ((1. - x) * (1 << 59 + csdp * (1. - x) * ((x + y - 1.) * ( << 60 + cif * (x - 1.) * (1. - y) / x + cig << 61 } << 62 << 63 const G4double d2wmax = D2W(xl, yl); << 64 } // namespace << 65 45 66 G4PionRadiativeDecayChannel::G4PionRadiativeDe << 46 G4PionRadiativeDecayChannel:: 67 << 47 G4PionRadiativeDecayChannel(const G4String& theParentName, 68 : G4VDecayChannel("Radiative Pion Decay", 1) << 48 G4double theBR) >> 49 : G4VDecayChannel("Radiative Pion Decay",1) 69 { 50 { 70 // set names for daughter particles 51 // set names for daughter particles 71 if (theParentName == "pi+") { 52 if (theParentName == "pi+") { 72 SetBR(theBR); 53 SetBR(theBR); 73 SetParent("pi+"); 54 SetParent("pi+"); 74 SetNumberOfDaughters(3); 55 SetNumberOfDaughters(3); 75 SetDaughter(0, "e+"); 56 SetDaughter(0, "e+"); 76 SetDaughter(1, "gamma"); 57 SetDaughter(1, "gamma"); 77 SetDaughter(2, "nu_e"); 58 SetDaughter(2, "nu_e"); 78 } << 59 } else if (theParentName == "pi-") { 79 else if (theParentName == "pi-") { << 80 SetBR(theBR); 60 SetBR(theBR); 81 SetParent("pi-"); 61 SetParent("pi-"); 82 SetNumberOfDaughters(3); 62 SetNumberOfDaughters(3); 83 SetDaughter(0, "e-"); 63 SetDaughter(0, "e-"); 84 SetDaughter(1, "gamma"); 64 SetDaughter(1, "gamma"); 85 SetDaughter(2, "anti_nu_e"); 65 SetDaughter(2, "anti_nu_e"); 86 } << 66 } else { 87 else { << 88 #ifdef G4VERBOSE 67 #ifdef G4VERBOSE 89 if (GetVerboseLevel() > 0) { << 68 if (GetVerboseLevel()>0) { 90 G4cout << "G4RadiativePionDecayChannel:: << 69 G4cout << "G4RadiativePionDecayChannel:: constructor :"; 91 G4cout << "Parent particle is not charge << 70 G4cout << " parent particle is not muon but "; 92 G4cout << theParentName << G4endl; 71 G4cout << theParentName << G4endl; 93 } 72 } 94 #endif 73 #endif 95 } 74 } >> 75 >> 76 beta = 3.6612e-03; >> 77 >> 78 cib = 1.16141e-03; >> 79 csdp = 3.45055e-02; >> 80 csdm = 5.14122e-03; >> 81 cif = 4.63543e-05; >> 82 cig = 1.78928e-05; >> 83 >> 84 xl = 2.*0.1*MeV/139.57*MeV; >> 85 yl = ((1.-xl) + std::sqrt((1-xl)*(1-xl)+4*beta*beta))/2.; >> 86 >> 87 xu = 1. - (yl - std::sqrt(yl*yl-4.*beta*beta))/2.; >> 88 yu = 1. + beta*beta; >> 89 >> 90 d2wmax = D2W(xl,yl); >> 91 96 } 92 } 97 93 98 G4PionRadiativeDecayChannel& << 94 G4PionRadiativeDecayChannel::~G4PionRadiativeDecayChannel() 99 G4PionRadiativeDecayChannel::operator=(const G << 100 { 95 { 101 if (this != &right) { << 102 kinematics_name = right.kinematics_name; << 103 verboseLevel = right.verboseLevel; << 104 rbranch = right.rbranch; << 105 << 106 // copy parent name << 107 parent_name = new G4String(*right.parent_n << 108 << 109 // clear daughters_name array << 110 ClearDaughtersName(); << 111 << 112 // recreate array << 113 numberOfDaughters = right.numberOfDaughter << 114 if (numberOfDaughters > 0) { << 115 if (daughters_name != nullptr) ClearDaug << 116 daughters_name = new G4String*[numberOfD << 117 // copy daughters name << 118 for (G4int index = 0; index < numberOfDa << 119 daughters_name[index] = new G4String(* << 120 } << 121 } << 122 } << 123 return *this; << 124 } 96 } 125 97 126 G4DecayProducts* G4PionRadiativeDecayChannel:: << 98 G4DecayProducts *G4PionRadiativeDecayChannel::DecayIt(G4double) 127 { 99 { >> 100 128 #ifdef G4VERBOSE 101 #ifdef G4VERBOSE 129 if (GetVerboseLevel() > 1) G4cout << "G4Pion << 102 if (GetVerboseLevel()>1) >> 103 G4cout << "G4PionRadiativeDecayChannel::DecayIt "; 130 #endif 104 #endif 131 105 132 CheckAndFillParent(); << 106 if (parent == 0) FillParent(); 133 CheckAndFillDaughters(); << 107 if (daughters == 0) FillDaughters(); 134 108 135 // parent mass 109 // parent mass 136 G4double parentmass = G4MT_parent->GetPDGMas << 110 G4double parentmass = parent->GetPDGMass(); 137 111 138 G4double EMPI = parentmass; 112 G4double EMPI = parentmass; 139 113 140 // daughters'mass << 114 //daughters'mass 141 const G4int N_DAUGHTER = 3; << 115 G4double daughtermass[3]; 142 G4double daughtermass[N_DAUGHTER]; << 116 G4double sumofdaughtermass = 0.0; 143 // G4double sumofdaughtermass = 0.0; << 117 for (G4int index=0; index<3; index++){ 144 for (G4int index = 0; index < N_DAUGHTER; ++ << 118 daughtermass[index] = daughters[index]->GetPDGMass(); 145 daughtermass[index] = G4MT_daughters[index << 119 sumofdaughtermass += daughtermass[index]; 146 // sumofdaughtermass += daughtermass[index << 147 } 120 } 148 121 149 G4double EMASS = daughtermass[0]; 122 G4double EMASS = daughtermass[0]; 150 123 151 // create parent G4DynamicParticle at rest << 124 //create parent G4DynamicParticle at rest 152 G4ThreeVector dummy; 125 G4ThreeVector dummy; 153 auto parentparticle = new G4DynamicParticle( << 126 G4DynamicParticle * parentparticle = 154 // create G4Decayproducts << 127 new G4DynamicParticle( parent, dummy, 0.0); 155 auto products = new G4DecayProducts(*parentp << 128 //create G4Decayproducts >> 129 G4DecayProducts *products = new G4DecayProducts(*parentparticle); 156 delete parentparticle; 130 delete parentparticle; 157 131 158 G4double x, y; << 132 G4double x, y, d2w; 159 133 160 const std::size_t MAX_LOOP = 1000; << 134 do { 161 135 162 for (std::size_t loop_counter1 = 0; loop_cou << 136 do { 163 for (std::size_t loop_counter2 = 0; loop_c << 137 164 x = xl + G4UniformRand() * (xu - xl); << 138 x = xl + G4UniformRand()*(xu-xl); 165 y = yl + G4UniformRand() * (yu - yl); << 139 y = yl + G4UniformRand()*(yu-yl); 166 if (x + y > 1.) break; << 140 167 } << 141 } while (x+y <= 1.); 168 G4double d2w = D2W(x, y); << 169 if (d2w > G4UniformRand() * d2wmax) break; << 170 } << 171 142 172 // Calculate the angle between positron and << 143 d2w = D2W(x,y); 173 // << 174 G4double cthetaGE = << 175 (y * (x - 2.) + 2. * (1. - x + beta * beta << 176 144 177 G4double G = x * EMPI / 2.; << 145 } while (d2w <= G4UniformRand()*d2wmax); 178 G4double E = y * EMPI / 2.; << 146 >> 147 //----------------------------------------------------------------------- >> 148 // >> 149 // Calculate the angle between positron and photon (cosine) >> 150 // >> 151 G4double cthetaGE = (y*(x-2.)+2.*(1.-x+beta*beta)) / >> 152 (x*std::sqrt(y*y-4.*beta*beta)); >> 153 >> 154 // >> 155 //----------------------------------------------------------------------- >> 156 // >> 157 G4double G = x * EMPI/2.; >> 158 G4double E = y * EMPI/2.; >> 159 // >> 160 //----------------------------------------------------------------------- >> 161 // 179 162 180 if (E < EMASS) E = EMASS; 163 if (E < EMASS) E = EMASS; 181 164 182 // calculate daughter momentum 165 // calculate daughter momentum 183 G4double daughtermomentum[2]; 166 G4double daughtermomentum[2]; 184 167 185 daughtermomentum[0] = std::sqrt(E * E - EMAS << 168 daughtermomentum[0] = std::sqrt(E*E - EMASS*EMASS); 186 169 187 G4double cthetaE = 2. * G4UniformRand() - 1. << 170 G4double cthetaE = 2.*G4UniformRand()-1.; 188 G4double sthetaE = std::sqrt(1. - cthetaE * << 171 G4double sthetaE = std::sqrt(1.-cthetaE*cthetaE); 189 172 190 G4double phiE = twopi * G4UniformRand() * ra << 173 G4double phiE = twopi*G4UniformRand()*rad; 191 G4double cphiE = std::cos(phiE); 174 G4double cphiE = std::cos(phiE); 192 G4double sphiE = std::sin(phiE); 175 G4double sphiE = std::sin(phiE); 193 176 194 // Coordinates of the decay positron << 177 //Coordinates of the decay positron 195 // << 178 196 G4double px = sthetaE * cphiE; << 179 G4double px = sthetaE*cphiE; 197 G4double py = sthetaE * sphiE; << 180 G4double py = sthetaE*sphiE; 198 G4double pz = cthetaE; 181 G4double pz = cthetaE; 199 182 200 G4ThreeVector direction0(px, py, pz); << 183 G4ThreeVector direction0(px,py,pz); 201 184 202 auto daughterparticle0 = << 185 G4DynamicParticle * daughterparticle0 203 new G4DynamicParticle(G4MT_daughters[0], d << 186 = new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0); 204 187 205 products->PushProducts(daughterparticle0); 188 products->PushProducts(daughterparticle0); 206 189 207 daughtermomentum[1] = G; 190 daughtermomentum[1] = G; 208 191 209 G4double sthetaGE = std::sqrt(1. - cthetaGE << 192 G4double sthetaGE = std::sqrt(1.-cthetaGE*cthetaGE); 210 193 211 G4double phiGE = twopi * G4UniformRand() * r << 194 G4double phiGE = twopi*G4UniformRand()*rad; 212 G4double cphiGE = std::cos(phiGE); 195 G4double cphiGE = std::cos(phiGE); 213 G4double sphiGE = std::sin(phiGE); 196 G4double sphiGE = std::sin(phiGE); 214 197 215 // Coordinates of the decay gamma with respe << 198 //Coordinates of the decay gamma with respect to the decay positron 216 // << 199 217 px = sthetaGE * cphiGE; << 200 px = sthetaGE*cphiGE; 218 py = sthetaGE * sphiGE; << 201 py = sthetaGE*sphiGE; 219 pz = cthetaGE; 202 pz = cthetaGE; 220 203 221 G4ThreeVector direction1(px, py, pz); << 204 G4ThreeVector direction1(px,py,pz); 222 205 223 direction1.rotateUz(direction0); 206 direction1.rotateUz(direction0); 224 207 225 auto daughterparticle1 = << 208 G4DynamicParticle * daughterparticle1 226 new G4DynamicParticle(G4MT_daughters[1], d << 209 = new G4DynamicParticle( daughters[1], daughtermomentum[1]*direction1); 227 210 228 products->PushProducts(daughterparticle1); 211 products->PushProducts(daughterparticle1); 229 212 230 // output message << 213 // output message 231 #ifdef G4VERBOSE 214 #ifdef G4VERBOSE 232 if (GetVerboseLevel() > 1) { << 215 if (GetVerboseLevel()>1) { 233 G4cout << "G4PionRadiativeDecayChannel::De << 216 G4cout << "G4PionRadiativeDecayChannel::DecayIt "; 234 G4cout << " create decay products in rest << 217 G4cout << " create decay products in rest frame " <<G4endl; 235 products->DumpInfo(); 218 products->DumpInfo(); 236 } 219 } 237 #endif 220 #endif 238 221 239 return products; 222 return products; >> 223 240 } 224 } 241 225