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 // ABLAXX statistical de-excitation model 26 // ABLAXX statistical de-excitation model 27 // Jose Luis Rodriguez, UDC (translation from << 27 // Pekka Kaitaniemi, HIP (translation) 28 // Pekka Kaitaniemi, HIP (initial translation << 28 // Christelle Schmidt, IPNL (fission code) 29 // Aleksandra Kelic, GSI (ABLA07 code) << 29 // Davide Mancusi, CEA (contact person INCL/ABLA) 30 // Davide Mancusi, CEA (contact person INCL) << 31 // Aatos Heikkinen, HIP (project coordination) 30 // Aatos Heikkinen, HIP (project coordination) 32 // 31 // >> 32 #define ABLAXX_IN_GEANT4_MODE 1 33 33 34 #include "globals.hh" 34 #include "globals.hh" 35 #include <cmath> << 35 36 #include <iostream> << 36 #ifdef ABLAXX_IN_GEANT4_MODE 37 37 38 #include "G4AblaInterface.hh" 38 #include "G4AblaInterface.hh" 39 #include "G4DoubleHyperDoubleNeutron.hh" << 40 #include "G4DoubleHyperH4.hh" << 41 #include "G4DynamicParticle.hh" << 42 #include "G4ExcitationHandler.hh" << 43 #include "G4HyperAlpha.hh" << 44 #include "G4HyperH4.hh" << 45 #include "G4HyperHe5.hh" << 46 #include "G4HyperTriton.hh" << 47 #include "G4IonTable.hh" << 48 #include "G4ParticleDefinition.hh" 39 #include "G4ParticleDefinition.hh" 49 #include "G4PhysicalConstants.hh" << 50 #include "G4PhysicsModelCatalog.hh" << 51 #include "G4ReactionProduct.hh" << 52 #include "G4ReactionProductVector.hh" 40 #include "G4ReactionProductVector.hh" >> 41 #include "G4ReactionProduct.hh" >> 42 #include "G4DynamicParticle.hh" >> 43 #include "G4IonTable.hh" 53 #include "G4SystemOfUnits.hh" 44 #include "G4SystemOfUnits.hh" >> 45 #include "G4PhysicalConstants.hh" >> 46 #include <iostream> >> 47 #include <cmath> 54 48 55 G4AblaInterface::G4AblaInterface(G4ExcitationH << 49 G4AblaInterface::G4AblaInterface() : 56 : G4VPreCompoundModel(ptr, "ABLAXX") << 50 G4VPreCompoundModel(NULL, "ABLA"), 57 , ablaResult(new G4VarNtp) << 51 ablaResult(new G4VarNtp), 58 , theABLAModel(new G4Abla(ablaResult)) << 52 volant(new G4Volant), 59 , eventNumber(0) << 53 theABLAModel(new G4Abla(volant, ablaResult)), 60 , secID(-1) << 54 eventNumber(0) 61 , isInitialised(false) << 55 { 62 { << 56 theABLAModel->initEvapora(); 63 secID = G4PhysicsModelCatalog::GetModelID( << 57 } 64 // G4cout << "### NEW PrecompoundModel " < << 58 65 if (!ptr) << 59 G4AblaInterface::~G4AblaInterface() { 66 SetExcitationHandler(new G4ExcitationH << 60 delete volant; 67 InitialiseModel(); << 61 delete ablaResult; 68 G4cout << G4endl << "G4AblaInterface::Init << 62 delete theABLAModel; 69 } << 63 } 70 << 64 71 G4AblaInterface::~G4AblaInterface() << 65 G4ReactionProductVector *G4AblaInterface::DeExcite(G4Fragment &aFragment) { 72 { << 66 volant->clear(); 73 applyYourselfResult.Clear(); << 67 ablaResult->clear(); 74 delete ablaResult; << 68 75 delete theABLAModel; << 69 const G4int ARem = aFragment.GetA_asInt(); 76 delete GetExcitationHandler(); << 70 const G4int ZRem = aFragment.GetZ_asInt(); 77 } << 71 const G4double nuclearMass = aFragment.GetGroundStateMass() / MeV; 78 << 72 const G4double eStarRem = aFragment.GetExcitationEnergy() / MeV; 79 void G4AblaInterface::BuildPhysicsTable(const << 73 const G4double jRem = aFragment.GetAngularMomentum().mag() / hbar_Planck; 80 << 74 const G4LorentzVector &pRem = aFragment.GetMomentum(); 81 void G4AblaInterface::InitialiseModel() << 75 const G4double eTotRem = pRem.e(); 82 { << 76 const G4double eKinRem = (eTotRem - pRem.invariantMass()) / MeV; 83 if (isInitialised) << 77 const G4double pxRem = pRem.x() / MeV; 84 return; << 78 const G4double pyRem = pRem.y() / MeV; 85 isInitialised = true; << 79 const G4double pzRem = pRem.z() / MeV; 86 theABLAModel->initEvapora(); << 80 87 theABLAModel->SetParameters(); << 81 eventNumber++; 88 GetExcitationHandler()->Initialise(); << 82 89 } << 83 theABLAModel->breakItUp(ARem, ZRem, 90 << 84 nuclearMass, 91 G4HadFinalState* G4AblaInterface::ApplyYoursel << 85 eStarRem, 92 { << 86 jRem, 93 // This method is adapted from G4PreCompo << 87 eKinRem, 94 // and it is used only by Binary Cascade ( << 88 pxRem, 95 // Abla for nuclear de-excitation. This me << 89 pyRem, 96 // for proton and neutron projectile with << 90 pzRem, 97 // creating a "compound" nucleus made by t << 91 eventNumber); 98 // projectile", before calling the DeExcit << 92 99 const G4ParticleDefinition* primary = theP << 93 G4ReactionProductVector *result = new G4ReactionProductVector; 100 if (primary != G4Neutron::Definition() && << 94 101 { << 95 for(int j = 0; j < ablaResult->ntrack; ++j) { // Copy ABLA result to the EventInfo 102 G4ExceptionDescription ed; << 96 G4ReactionProduct *product = toG4Particle(ablaResult->avv[j], 103 ed << "G4AblaModel is used for "; << 97 ablaResult->zvv[j], 104 if (primary) << 98 ablaResult->enerj[j], 105 ed << primary->GetParticleName(); << 99 ablaResult->plab[j]*std::sin(ablaResult->tetlab[j]*pi/180.0)*std::cos(ablaResult->philab[j]*pi/180.0), 106 G4Exception("G4AblaInterface::ApplyYou << 100 ablaResult->plab[j]*std::sin(ablaResult->tetlab[j]*pi/180.0)*std::sin(ablaResult->philab[j]*pi/180.0), 107 return nullptr; << 101 ablaResult->plab[j]*std::cos(ablaResult->tetlab[j]*pi/180.0)); 108 } << 102 if(product) 109 << 103 result->push_back(product); 110 G4int Zp = 0; << 104 } 111 G4int Ap = 1; << 105 return result; 112 if (primary == G4Proton::Definition()) << 106 } 113 Zp = 1; << 107 114 G4double timePrimary = thePrimary.GetGloba << 108 G4ParticleDefinition *G4AblaInterface::toG4ParticleDefinition(G4int A, G4int Z) const { 115 G4int A = theNucleus.GetA_asInt(); << 109 if (A == 1 && Z == 1) return G4Proton::Proton(); 116 G4int Z = theNucleus.GetZ_asInt(); << 110 else if(A == 1 && Z == 0) return G4Neutron::Neutron(); 117 G4LorentzVector p = thePrimary.Get4Momentu << 111 else if(A == 0 && Z == 1) return G4PionPlus::PionPlus(); 118 G4double mass = G4NucleiProperties::GetNuc << 112 else if(A == 0 && Z == -1) return G4PionMinus::PionMinus(); 119 p += G4LorentzVector(0.0, 0.0, 0.0, mass); << 113 else if(A == 0 && Z == 0) return G4PionZero::PionZero(); 120 << 114 else if(A == 2 && Z == 1) return G4Deuteron::Deuteron(); 121 G4Fragment anInitialState(A + Ap, Z + Zp, << 115 else if(A == 3 && Z == 1) return G4Triton::Triton(); 122 anInitialState.SetNumberOfExcitedParticle( << 116 else if(A == 3 && Z == 2) return G4He3::He3(); 123 anInitialState.SetNumberOfHoles(1, Zp); << 117 else if(A == 4 && Z == 2) return G4Alpha::Alpha(); 124 anInitialState.SetCreationTime(thePrimary. << 118 else if(A > 0 && Z > 0 && A >= Z) { // Returns ground state ion definition 125 anInitialState.SetCreatorModelID(secID); << 119 return G4IonTable::GetIonTable()->GetIon(Z, A, 0); 126 << 120 } else { // Error, unrecognized particle 127 G4ReactionProductVector* deExciteResult = << 121 G4cout << "Can't convert particle with A=" << A << ", Z=" << Z << " to G4ParticleDefinition, trouble ahead" << G4endl; 128 << 122 return 0; 129 applyYourselfResult.Clear(); << 123 } 130 applyYourselfResult.SetStatusChange(stopAn << 124 } 131 for (auto const& prod : *deExciteResult) << 125 132 { << 126 G4ReactionProduct *G4AblaInterface::toG4Particle(G4int A, G4int Z, 133 G4DynamicParticle* aNewDP = << 127 G4double kinE, 134 new G4DynamicParticle(prod->GetDef << 128 G4double px, 135 G4HadSecondary aNew = G4HadSecondary(a << 129 G4double py, G4double pz) const { 136 G4double time = std::max(prod->GetForm << 130 G4ParticleDefinition *def = toG4ParticleDefinition(A, Z); 137 aNew.SetTime(timePrimary + time); << 131 if(def == 0) { // Check if we have a valid particle definition 138 aNew.SetCreatorModelID(prod->GetCreato << 132 return 0; 139 delete prod; << 133 } 140 applyYourselfResult.AddSecondary(aNew) << 134 const double energy = kinE * MeV; 141 } << 135 const G4ThreeVector momentum(px, py, pz); 142 delete deExciteResult; << 136 const G4ThreeVector momentumDirection = momentum.unit(); 143 return &applyYourselfResult; << 137 G4DynamicParticle p(def, momentumDirection, energy); 144 } << 138 G4ReactionProduct *r = new G4ReactionProduct(def); 145 << 139 (*r) = p; 146 G4ReactionProductVector* G4AblaInterface::DeEx << 140 return r; 147 { << 141 } 148 if (!isInitialised) << 142 149 InitialiseModel(); << 143 void G4AblaInterface::ModelDescription(std::ostream& outFile) const { 150 << 144 outFile << "ABLA V3 does not provide an implementation of the ApplyYourself method!\n\n"; 151 ablaResult->clear(); << 145 } 152 << 146 153 const G4int ARem = aFragment.GetA_asInt(); << 147 void G4AblaInterface::DeExciteModelDescription(std::ostream& outFile) const { 154 const G4int ZRem = aFragment.GetZ_asInt(); << 148 outFile << "ABLA V3 is a statistical model for nuclear de-excitation. It simulates\n" 155 const G4int SRem = -aFragment.GetNumberOfL << 149 << "evaporation of neutrons, protons and alpha particles, as well as fission\n" 156 const G4double eStarRem = aFragment.GetExc << 150 << "where applicable. The code included in Geant4 is a C++ translation of the\n" 157 const G4double jRem = aFragment.GetAngular << 151 << "original Fortran code. More details about the physics are available in the\n" 158 const G4LorentzVector& pRem = aFragment.Ge << 152 << "the Geant4 Physics Reference Manual and in the reference articles.\n\n" 159 const G4double pxRem = pRem.x() / MeV; << 153 << "References: A.R. Junghans et al., Nucl. Phys. A629 (1998) 635;\n" 160 const G4double pyRem = pRem.y() / MeV; << 154 << " J. Benlliure et al., Nucl. Phys. A628 (1998) 458.\n\n"; } 161 const G4double pzRem = pRem.z() / MeV; << 162 << 163 ++eventNumber; << 164 << 165 theABLAModel->DeexcitationAblaxx(ARem, ZRe << 166 << 167 G4ReactionProductVector* result = new G4Re << 168 << 169 for (G4int j = 0; j < ablaResult->ntrack; << 170 { // Copy ABLA result to the EventInfo << 171 G4ReactionProduct* product = toG4Parti << 172 << 173 << 174 << 175 << 176 << 177 << 178 if (product) << 179 { << 180 product->SetCreatorModelID(secID); << 181 result->push_back(product); << 182 } << 183 } << 184 return result; << 185 } << 186 << 187 G4ParticleDefinition* G4AblaInterface::toG4Par << 188 { << 189 if (A == 1 && Z == 1 && S == 0) << 190 return G4Proton::Proton(); << 191 else if (A == 1 && Z == 0 && S == 0) << 192 return G4Neutron::Neutron(); << 193 else if (A == 1 && Z == 0 && S == -1) << 194 return G4Lambda::Lambda(); << 195 else if (A == -1 && Z == 1 && S == 0) << 196 return G4PionPlus::PionPlus(); << 197 else if (A == -1 && Z == -1 && S == 0) << 198 return G4PionMinus::PionMinus(); << 199 else if (A == -1 && Z == 0 && S == 0) << 200 return G4PionZero::PionZero(); << 201 else if (A == 0 && Z == 0 && S == 0) << 202 return G4Gamma::Gamma(); << 203 else if (A == 2 && Z == 1 && S == 0) << 204 return G4Deuteron::Deuteron(); << 205 else if (A == 3 && Z == 1 && S == 0) << 206 return G4Triton::Triton(); << 207 else if (A == 3 && Z == 2 && S == 0) << 208 return G4He3::He3(); << 209 else if (A == 3 && Z == 1 && S == -1) << 210 return G4HyperTriton::Definition(); << 211 else if (A == 4 && Z == 2 && S == 0) << 212 return G4Alpha::Alpha(); << 213 else if (A == 4 && Z == 1 && S == -1) << 214 return G4HyperH4::Definition(); << 215 else if (A == 4 && Z == 2 && S == -1) << 216 return G4HyperAlpha::Definition(); << 217 else if (A == 4 && Z == 1 && S == -2) << 218 return G4DoubleHyperH4::Definition(); << 219 else if (A == 4 && Z == 0 && S == -2) << 220 return G4DoubleHyperDoubleNeutron::Def << 221 else if (A == 5 && Z == 2 && S == -1) << 222 return G4HyperHe5::Definition(); << 223 else if (A > 0 && Z > 0 && A > Z) << 224 { // Returns ground state ion definition. << 225 auto ionfromtable = G4IonTable::GetIon << 226 if (ionfromtable) << 227 return ionfromtable; << 228 else << 229 { << 230 G4cout << "Can't convert particle << 231 << " to G4ParticleDefinitio << 232 return 0; << 233 } << 234 } << 235 else << 236 { // Error, unrecognized particle << 237 G4cout << "Can't convert particle with << 238 << " to G4ParticleDefinition, t << 239 return 0; << 240 } << 241 } << 242 << 243 G4ReactionProduct* << 244 G4AblaInterface::toG4Particle(G4int A, G4i << 245 { << 246 G4ParticleDefinition* def = toG4ParticleDe << 247 if (def == 0) << 248 { // Check if we have a valid particle def << 249 return 0; << 250 } << 251 << 252 const G4double energy = kinE * MeV; << 253 const G4ThreeVector momentum(px, py, pz); << 254 const G4ThreeVector momentumDirection = mo << 255 G4DynamicParticle p(def, momentumDirection << 256 G4ReactionProduct* r = new G4ReactionProdu << 257 (*r) = p; << 258 return r; << 259 } << 260 << 261 void G4AblaInterface::ModelDescription(std::os << 262 { << 263 outFile << "ABLA++ does not provide an imp << 264 "method!\n\n"; << 265 } << 266 155 267 void G4AblaInterface::DeExciteModelDescription << 156 #endif // ABLAXX_IN_GEANT4_MODE 268 { << 269 outFile << "ABLA++ is a statistical model << 270 << "the gamma emission and the eva << 271 << "particles and IMFs, as well as << 272 << "included in Geant4 is a C++ tr << 273 << "code ABLA07. Although the mode << 274 << "hypernuclei by including the e << 275 << "More details about the physics << 276 << "Physics Reference Manual and i << 277 << "References:\n" << 278 << "(1) A. Kelic, M. V. Ricciardi, << 279 << "ICTP-IAEA Advanced Workshop on << 280 << "ICTP Trieste, Italy, 4–8 Feb << 281 "Leray, Y. Yariv, A. Mengoni, A << 282 "INDC(NDS)-530, Vienna, 2008), << 283 << "(2) J.L. Rodriguez-Sanchez, J. << 284 << "(3) J.L. Rodriguez-Sanchez et << 285 << "(4) J.L. Rodriguez-Sanchez et << 286 } << 287 157