00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 #include <assert.h>
00054 #include <sstream>
00055 #include <map>
00056
00057 #include "G4InuclNuclei.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "G4Fragment.hh"
00060 #include "G4HadronicException.hh"
00061 #include "G4InuclSpecialFunctions.hh"
00062 #include "G4Ions.hh"
00063 #include "G4IonTable.hh"
00064 #include "G4NucleiProperties.hh"
00065 #include "G4Nucleon.hh"
00066 #include "G4ParticleDefinition.hh"
00067 #include "G4ParticleTable.hh"
00068 #include "G4V3DNucleus.hh"
00069
00070 using namespace G4InuclSpecialFunctions;
00071
00072
00073
00074
00075 G4InuclNuclei::G4InuclNuclei(const G4Fragment& aFragment,
00076 G4InuclParticle::Model model)
00077 : G4InuclParticle() {
00078 copy(aFragment, model);
00079 }
00080
00081 void G4InuclNuclei::copy(const G4Fragment& aFragment, Model model) {
00082 fill(aFragment.GetMomentum()/GeV, aFragment.GetA_asInt(),
00083 aFragment.GetZ_asInt(), aFragment.GetExcitationEnergy(), model);
00084
00085
00086 theExitonConfiguration.protonQuasiParticles = aFragment.GetNumberOfCharged();
00087
00088 theExitonConfiguration.neutronQuasiParticles =
00089 aFragment.GetNumberOfParticles() - aFragment.GetNumberOfCharged();
00090
00091 theExitonConfiguration.protonHoles = aFragment.GetNumberOfChargedHoles();
00092
00093 theExitonConfiguration.neutronHoles =
00094 aFragment.GetNumberOfHoles() - theExitonConfiguration.protonHoles;
00095 }
00096
00097
00098
00099 G4Fragment G4InuclNuclei::makeG4Fragment() const {
00100 G4Fragment frag(getA(), getZ(), getMomentum()*GeV);
00101
00102
00103 frag.SetNumberOfHoles(theExitonConfiguration.protonHoles
00104 + theExitonConfiguration.neutronHoles,
00105 theExitonConfiguration.protonHoles);
00106
00107 frag.SetNumberOfExcitedParticle(theExitonConfiguration.protonQuasiParticles
00108 + theExitonConfiguration.neutronQuasiParticles,
00109 theExitonConfiguration.protonQuasiParticles);
00110
00111 return frag;
00112 }
00113
00114 G4InuclNuclei::operator G4Fragment() const {
00115 return makeG4Fragment();
00116 }
00117
00118
00119
00120
00121 G4InuclNuclei::G4InuclNuclei(G4V3DNucleus* a3DNucleus,
00122 G4InuclParticle::Model model)
00123 : G4InuclParticle() {
00124 copy(a3DNucleus, model);
00125 }
00126
00127 void G4InuclNuclei::copy(G4V3DNucleus* a3DNucleus, Model model) {
00128 if (!a3DNucleus) return;
00129
00130 fill(0., a3DNucleus->GetMassNumber(), a3DNucleus->GetCharge(), 0., model);
00131
00132
00133 if (a3DNucleus->StartLoop()) {
00134 G4Nucleon* nucl = 0;
00135 while ((nucl = a3DNucleus->GetNextNucleon())) {
00136 if (nucl->AreYouHit()) {
00137 if (nucl->GetParticleType() == G4Proton::Definition())
00138 theExitonConfiguration.protonHoles++;
00139
00140 if (nucl->GetParticleType() == G4Neutron::Definition())
00141 theExitonConfiguration.neutronHoles++;
00142 }
00143 }
00144 }
00145 }
00146
00147
00148
00149
00150 void G4InuclNuclei::fill(const G4LorentzVector& mom, G4int a, G4int z,
00151 G4double exc, G4InuclParticle::Model model) {
00152 setDefinition(makeDefinition(a,z));
00153 setMomentum(mom);
00154 setExitationEnergy(exc);
00155 clearExitonConfiguration();
00156 setModel(model);
00157 }
00158
00159 void G4InuclNuclei::fill(G4double ekin, G4int a, G4int z, G4double exc,
00160 G4InuclParticle::Model model) {
00161 setDefinition(makeDefinition(a,z));
00162 setKineticEnergy(ekin);
00163 setExitationEnergy(exc);
00164 clearExitonConfiguration();
00165 setModel(model);
00166 }
00167
00168 void G4InuclNuclei::clear() {
00169 setDefinition(0);
00170 clearExitonConfiguration();
00171 setModel(G4InuclParticle::DefaultModel);
00172 }
00173
00174
00175
00176
00177 void G4InuclNuclei::setExitationEnergy(G4double e) {
00178 G4double ekin = getKineticEnergy();
00179
00180 G4double emass = getNucleiMass() + e*MeV/GeV;
00181
00182
00183 G4double ekin_new = std::sqrt(emass*emass + ekin*(2.*getMass()+ekin)) - emass;
00184
00185 setMass(emass);
00186 setKineticEnergy(ekin_new);
00187 }
00188
00189
00190
00191
00192
00193
00194
00195 G4ParticleDefinition* G4InuclNuclei::makeDefinition(G4int a, G4int z) {
00196
00197 if (0 == a && 0 == z) return 0;
00198
00199 G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
00200 G4ParticleDefinition *pd = pTable->GetIon(z, a, 0.);
00201
00202
00203 if (!pd) pd = makeNuclearFragment(a,z);
00204
00205 return pd;
00206 }
00207
00208
00209
00210
00211
00212
00213 G4ParticleDefinition*
00214 G4InuclNuclei::makeNuclearFragment(G4int a, G4int z) {
00215 if (a<=0 || z<0 || a<z) {
00216 G4cerr << " >>> G4InuclNuclei::makeNuclearFragment() called with"
00217 << " impossible arguments A=" << a << " Z=" << z << G4endl;
00218 throw G4HadronicException(__FILE__, __LINE__,
00219 "G4InuclNuclei impossible A/Z arguments");
00220 }
00221
00222 G4int code = G4IonTable::GetNucleusEncoding(z, a);
00223
00224
00225
00226
00227
00228
00229 static std::map<G4int, G4ParticleDefinition*> fragmentList;
00230 if (fragmentList.find(code) != fragmentList.end()) return fragmentList[code];
00231
00232
00233 std::stringstream zstr, astr;
00234 zstr << z;
00235 astr << a;
00236
00237 G4String name = "Z" + zstr.str() + "A" + astr.str();
00238
00239 G4double mass = getNucleiMass(a,z) *GeV/MeV;
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249 G4Ions* fragPD = new G4Ions(name, mass, 0., z*eplus,
00250 0, +1, 0,
00251 0, 0, 0,
00252 "nucleus", 0, a, code,
00253 true, 0., 0,
00254 true, "generic", 0, 0.);
00255 fragPD->SetAntiPDGEncoding(0);
00256
00257 return (fragmentList[code] = fragPD);
00258 }
00259
00260 G4double G4InuclNuclei::getNucleiMass(G4int a, G4int z, G4double exc) {
00261
00262 G4double mass = G4NucleiProperties::GetNuclearMass(a,z) + exc;
00263
00264 return mass*MeV/GeV;
00265 }
00266
00267
00268 G4InuclNuclei& G4InuclNuclei::operator=(const G4InuclNuclei& right) {
00269 theExitonConfiguration = right.theExitonConfiguration;
00270 G4InuclParticle::operator=(right);
00271 return *this;
00272 }
00273
00274
00275
00276 void G4InuclNuclei::print(std::ostream& os) const {
00277 G4InuclParticle::print(os);
00278 os << G4endl << " Nucleus: " << getDefinition()->GetParticleName()
00279 << " A " << getA() << " Z " << getZ() << " mass " << getMass()
00280 << " Eex (MeV) " << getExitationEnergy();
00281
00282 if (!theExitonConfiguration.empty())
00283 os << G4endl << " " << theExitonConfiguration;
00284 }