Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 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 // 26 // 27 // 28 // P. Arce, June-2014 Conversion neutron_hp to 29 // 30 #include "G4ParticleHPNBodyPhaseSpace.hh" 31 32 #include "G4Alpha.hh" 33 #include "G4Deuteron.hh" 34 #include "G4Electron.hh" 35 #include "G4Gamma.hh" 36 #include "G4He3.hh" 37 #include "G4Neutron.hh" 38 #include "G4PhysicalConstants.hh" 39 #include "G4Positron.hh" 40 #include "G4Proton.hh" 41 #include "G4ThreeVector.hh" 42 #include "G4Triton.hh" 43 #include "Randomize.hh" 44 45 G4ReactionProduct* G4ParticleHPNBodyPhaseSpace 46 47 { 48 auto result = new G4ReactionProduct; 49 auto Z = static_cast<G4int>(massCode / 1000) 50 auto A = static_cast<G4int>(massCode - 1000 51 52 if (massCode == 0) { 53 result->SetDefinition(G4Gamma::Gamma()); 54 } 55 else if (A == 0) { 56 result->SetDefinition(G4Electron::Electron 57 if (Z == 1) result->SetDefinition(G4Positr 58 } 59 else if (A == 1) { 60 result->SetDefinition(G4Neutron::Neutron() 61 if (Z == 1) result->SetDefinition(G4Proton 62 } 63 else if (A == 2) { 64 result->SetDefinition(G4Deuteron::Deuteron 65 } 66 else if (A == 3) { 67 result->SetDefinition(G4Triton::Triton()); 68 if (Z == 2) result->SetDefinition(G4He3::H 69 } 70 else if (A == 4) { 71 result->SetDefinition(G4Alpha::Alpha()); 72 if (Z != 2) throw G4HadronicException(__FI 73 } 74 else { 75 throw G4HadronicException(__FILE__, __LINE 76 "G4ParticleHPNBo 77 } 78 79 // Get the energy from phase-space distribut 80 // in CMS 81 // P = Cn*std::sqrt(E')*(Emax-E')**(3*n/2-4) 82 G4double maxE = GetEmax(anEnergy, result->Ge 83 if (maxE <= 0) { 84 maxE = 1. * CLHEP::eV; 85 } 86 G4double energy = 0.; 87 G4double max(0); 88 if (theTotalCount <= 3) { 89 max = maxE / 2.; 90 } 91 else if (theTotalCount == 4) { 92 max = maxE / 5.; 93 } 94 else if (theTotalCount == 5) { 95 max = maxE / 8.; 96 } 97 else { 98 throw G4HadronicException( 99 __FILE__, __LINE__, 100 "NeutronHP Phase-space distribution cann 101 } 102 G4double testit; 103 G4double rand0 = Prob(max, maxE, theTotalCou 104 G4double rand; 105 106 G4int icounter = 0; 107 G4int icounter_max = 1024; 108 do { 109 icounter++; 110 if (icounter > icounter_max) { 111 G4cout << "Loop-counter exceeded the thr 112 << __FILE__ << "." << G4endl; 113 break; 114 } 115 rand = rand0 * G4UniformRand(); 116 energy = maxE * G4UniformRand(); 117 testit = Prob(energy, maxE, theTotalCount) 118 } while (rand > testit); // Loop checking, 119 result->SetKineticEnergy(energy); 120 121 // now do random direction 122 G4double cosTh = 2. * G4UniformRand() - 1.; 123 G4double phi = twopi * G4UniformRand(); 124 G4double theta = std::acos(cosTh); 125 G4double sinth = std::sin(theta); 126 G4double mtot = result->GetTotalMomentum(); 127 G4ThreeVector tempVector(mtot * sinth * std: 128 mtot * std::cos(the 129 result->SetMomentum(tempVector); 130 G4ReactionProduct aCMS = *GetTarget() + *Get 131 result->Lorentz(*result, -1. * aCMS); 132 return result; 133 } 134