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 #include <algorithm>
00040 #include <vector>
00041
00042 #include "G4GeneratorPrecompoundInterface.hh"
00043 #include "G4PhysicalConstants.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4DynamicParticleVector.hh"
00046 #include "G4KineticTrackVector.hh"
00047 #include "G4Proton.hh"
00048 #include "G4Neutron.hh"
00049 #include "G4V3DNucleus.hh"
00050 #include "G4Nucleon.hh"
00051 #include "G4FragmentVector.hh"
00052 #include "G4ReactionProduct.hh"
00053 #include "G4ReactionProductVector.hh"
00054 #include "G4PreCompoundModel.hh"
00055 #include "G4ExcitationHandler.hh"
00056 #include "G4DecayKineticTracks.hh"
00057 #include "G4HadronicInteractionRegistry.hh"
00058
00059 G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface(G4VPreCompoundModel* preModel)
00060 : CaptureThreshold(10*MeV)
00061 {
00062 proton = G4Proton::Proton();
00063 neutron = G4Neutron::Neutron();
00064 if(preModel) { SetDeExcitation(preModel); }
00065 else {
00066 G4HadronicInteraction* hadi =
00067 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
00068 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(hadi);
00069 if(!pre) { pre = new G4PreCompoundModel(); }
00070 SetDeExcitation(pre);
00071 }
00072 }
00073
00074 G4GeneratorPrecompoundInterface::~G4GeneratorPrecompoundInterface()
00075 {
00076 }
00077
00078
00079
00080 #define exactExcitationEnergy
00081
00082
00083
00084 G4ReactionProductVector* G4GeneratorPrecompoundInterface::
00085 Propagate(G4KineticTrackVector* theSecondaries, G4V3DNucleus* theNucleus)
00086 {
00087 G4ReactionProductVector * theTotalResult = new G4ReactionProductVector;
00088
00089
00090 G4DecayKineticTracks decay(theSecondaries);
00091
00092
00093 G4int anA=theNucleus->GetMassNumber();
00094 G4int aZ=theNucleus->GetCharge();
00095 G4int numberOfEx = 0;
00096 G4int numberOfCh = 0;
00097 G4int numberOfHoles = 0;
00098 G4double exEnergy = 0.0;
00099 G4double R = theNucleus->GetNuclearRadius();
00100 G4ThreeVector exciton3Momentum(0.,0.,0.);
00101 G4ThreeVector captured3Momentum(0.,0.,0.);
00102 G4ThreeVector wounded3Momentum(0.,0.,0.);
00103
00104
00105 #ifdef exactExcitationEnergy
00106 G4LorentzVector secondary4Momemtum(0,0,0,0);
00107 #endif
00108 G4KineticTrackVector::iterator iter;
00109 for(iter=theSecondaries->begin(); iter !=theSecondaries->end(); ++iter)
00110 {
00111 G4ParticleDefinition* part = (*iter)->GetDefinition();
00112 G4double e = (*iter)->Get4Momentum().e();
00113 G4double mass = (*iter)->Get4Momentum().mag();
00114 G4ThreeVector mom = (*iter)->Get4Momentum().vect();
00115 if((part != proton && part != neutron) ||
00116 (e > mass + CaptureThreshold) ||
00117 ((*iter)->GetPosition().mag() > R)) {
00118 G4ReactionProduct * theNew = new G4ReactionProduct(part);
00119 theNew->SetMomentum(mom);
00120 theNew->SetTotalEnergy(e);
00121 theTotalResult->push_back(theNew);
00122 #ifdef exactExcitationEnergy
00123 secondary4Momemtum += (*iter)->Get4Momentum();
00124 #endif
00125 } else {
00126
00127
00128 ++anA;
00129 ++numberOfEx;
00130 G4int Z = G4int(part->GetPDGCharge()/eplus + 0.1);
00131 aZ += Z;
00132 numberOfCh += Z;
00133 captured3Momentum += mom;
00134 exEnergy += (e - mass);
00135 }
00136 delete (*iter);
00137 }
00138 delete theSecondaries;
00139
00140
00141 G4Nucleon * theCurrentNucleon =
00142 theNucleus->StartLoop() ? theNucleus->GetNextNucleon() : 0;
00143 while(theCurrentNucleon) {
00144 if(theCurrentNucleon->AreYouHit()) {
00145 ++numberOfHoles;
00146 ++numberOfEx;
00147 --anA;
00148 aZ -= G4int(theCurrentNucleon->GetDefinition()->GetPDGCharge()/eplus + 0.1);
00149 wounded3Momentum += theCurrentNucleon->Get4Momentum().vect();
00150
00151 exEnergy += theCurrentNucleon->GetBindingEnergy();
00152 }
00153 theCurrentNucleon = theNucleus->GetNextNucleon();
00154 }
00155 exciton3Momentum = captured3Momentum - wounded3Momentum;
00156
00157 if(anA>0 && aZ>0) {
00158 G4double fMass = G4NucleiProperties::GetNuclearMass(anA, aZ);
00159 #ifdef exactExcitationEnergy
00160
00161 const G4HadProjectile * primary = GetPrimaryProjectile();
00162 G4double Einitial= primary->Get4Momentum().e()
00163 + G4NucleiProperties::GetNuclearMass(theNucleus->GetMassNumber(),theNucleus->GetCharge());
00164 G4double Efinal = fMass + secondary4Momemtum.e();
00165 if ( (Einitial - Efinal) > 0 ) {
00166
00167
00168 exEnergy=Einitial - Efinal;
00169 }
00170 else {
00171
00172
00173 exEnergy=0.;
00174 }
00175 #endif
00176 fMass += exEnergy;
00177 G4ThreeVector balance=primary->Get4Momentum().vect() - secondary4Momemtum.vect() - exciton3Momentum;
00178 #ifdef G4GPI_debug_excitation
00179 G4cout << "momentum balance init/final " << balance << " value " << balance.mag() << G4endl
00180 << "primary / secondaries "<< primary->Get4Momentum() << " / "
00181 << secondary4Momemtum << " captured/wounded: " << captured3Momentum << " / " << wounded3Momentum
00182 << " exciton " << exciton3Momentum << G4endl
00183 << secondary4Momemtum.vect() + exciton3Momentum << G4endl;
00184 #endif
00185 #ifdef exactExcitationEnergy
00186 G4LorentzVector exciton4Momentum(exciton3Momentum, fMass);
00187 #else
00188 G4LorentzVector exciton4Momentum(exciton3Momentum,
00189 std::sqrt(exciton3Momentum.mag2() + fMass*fMass));
00190 #endif
00191 if ( exEnergy > 0.0 ) {
00192 G4Fragment anInitialState(anA, aZ, exciton4Momentum);
00193 anInitialState.SetNumberOfParticles(numberOfEx-numberOfHoles);
00194 anInitialState.SetNumberOfCharged(numberOfCh);
00195 anInitialState.SetNumberOfHoles(numberOfHoles);
00196
00197 G4ReactionProductVector * aPrecoResult = theDeExcitation->DeExcite(anInitialState);
00198
00199 theTotalResult->insert(theTotalResult->end(),aPrecoResult->begin(),aPrecoResult->end() );
00200 delete aPrecoResult;
00201
00202 } else {
00203
00204 G4ParticleDefinition* theKindOfFragment = 0;
00205 if (anA == 1 && aZ == 0) {
00206 theKindOfFragment = G4Neutron::NeutronDefinition();
00207 } else if (anA == 1 && aZ == 1) {
00208 theKindOfFragment = G4Proton::ProtonDefinition();
00209 } else if (anA == 2 && aZ == 1) {
00210 theKindOfFragment = G4Deuteron::DeuteronDefinition();
00211 } else if (anA == 3 && aZ == 1) {
00212 theKindOfFragment = G4Triton::TritonDefinition();
00213 } else if (anA == 3 && aZ == 2) {
00214 theKindOfFragment = G4He3::He3Definition();
00215 } else if (anA == 4 && aZ == 2) {
00216 theKindOfFragment = G4Alpha::AlphaDefinition();;
00217 } else {
00218 theKindOfFragment =
00219 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(aZ,anA,0.0);
00220 }
00221 if (theKindOfFragment != 0) {
00222 G4ReactionProduct * theNew = new G4ReactionProduct(theKindOfFragment);
00223 theNew->SetMomentum(G4ThreeVector(0.,0.,0.));
00224 theNew->SetTotalEnergy(fMass);
00225
00226 theTotalResult->push_back(theNew);
00227 }
00228 }
00229 }
00230
00231 return theTotalResult;
00232 }
00233
00234 G4HadFinalState* G4GeneratorPrecompoundInterface::
00235 ApplyYourself(const G4HadProjectile &, G4Nucleus & )
00236 {
00237 G4cout << "G4GeneratorPrecompoundInterface: ApplyYourself interface called stand-allone."
00238 << G4endl;
00239 G4cout << "This class is only a mediator between generator and precompound"<<G4endl;
00240 G4cout << "Please remove from your physics list."<<G4endl;
00241 throw G4HadronicException(__FILE__, __LINE__, "SEVERE: G4GeneratorPrecompoundInterface model interface called stand-allone.");
00242 return new G4HadFinalState;
00243 }
00244 void G4GeneratorPrecompoundInterface::PropagateModelDescription(std::ostream& outFile) const
00245 {
00246 outFile << "G4GeneratorPrecompoundInterface interfaces a high\n"
00247 << "energy model through the wounded nucleus to precompound de-excition.\n"
00248 << "Low energy protons and neutron present among secondaries produced by \n"
00249 << "the high energy generator and within the nucleus are captured. The wounded\n"
00250 << "nucleus and the captured particles form an excited nuclear fragment. This\n"
00251 << "fragment is passed to the Geant4 pre-compound model for de-excitation.\n"
00252 << "Nuclear de-excitation:\n";
00253
00254
00255 }