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 #include <algorithm>
00041
00042 #include "G4LowEIonFragmentation.hh"
00043 #include "G4PhysicalConstants.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4Fancy3DNucleus.hh"
00046 #include "G4Proton.hh"
00047 #include "G4NucleiProperties.hh"
00048
00049 G4int G4LowEIonFragmentation::hits = 0;
00050 G4int G4LowEIonFragmentation::totalTries = 0;
00051 G4double G4LowEIonFragmentation::area = 0;
00052
00053 G4LowEIonFragmentation::G4LowEIonFragmentation(G4ExcitationHandler * const value)
00054 {
00055 theHandler = value;
00056 theModel = new G4PreCompoundModel(theHandler);
00057 proton = G4Proton::Proton();
00058 }
00059
00060 G4LowEIonFragmentation::G4LowEIonFragmentation()
00061 {
00062 theHandler = new G4ExcitationHandler;
00063 theModel = new G4PreCompoundModel(theHandler);
00064 proton = G4Proton::Proton();
00065 }
00066
00067 G4LowEIonFragmentation::~G4LowEIonFragmentation()
00068 {
00069 delete theModel;
00070 }
00071
00072 G4HadFinalState * G4LowEIonFragmentation::
00073 ApplyYourself(const G4HadProjectile & thePrimary, G4Nucleus & theNucleus)
00074 {
00075 area = 0;
00076
00077 theResult.Clear();
00078 theResult.SetStatusChange( stopAndKill );
00079 theResult.SetEnergyChange( 0.0 );
00080
00081
00082 G4int aTargetA = theNucleus.GetA_asInt();
00083 G4int aTargetZ = theNucleus.GetZ_asInt();
00084
00085
00086 G4int aProjectileA = thePrimary.GetDefinition()->GetBaryonNumber();
00087 G4int aProjectileZ = G4lrint(thePrimary.GetDefinition()->GetPDGCharge()/eplus);
00088
00089
00090
00091 G4Fancy3DNucleus aPrim;
00092 aPrim.Init(aProjectileA, aProjectileZ);
00093 G4double projectileOuterRadius = aPrim.GetOuterRadius();
00094
00095 G4Fancy3DNucleus aTarg;
00096 aTarg.Init(aTargetA, aTargetZ);
00097 G4double targetOuterRadius = aTarg.GetOuterRadius();
00098
00099
00100 G4int particlesFromProjectile = 0;
00101 G4int chargedFromProjectile = 0;
00102 G4double impactParameter = 0;
00103 G4double x,y;
00104 G4Nucleon * pNucleon;
00105
00106
00107 while(0==particlesFromProjectile)
00108 {
00109 do
00110 {
00111 x = 2*G4UniformRand() - 1;
00112 y = 2*G4UniformRand() - 1;
00113 }
00114 while(x*x + y*y > 1);
00115 impactParameter = std::sqrt(x*x+y*y)*(targetOuterRadius+projectileOuterRadius);
00116 ++totalTries;
00117 area = pi*(targetOuterRadius+projectileOuterRadius)*
00118 (targetOuterRadius+projectileOuterRadius);
00119 G4double projectileHorizon = impactParameter-targetOuterRadius;
00120
00121
00122 G4double empirical = G4UniformRand();
00123 if(projectileHorizon > empirical*projectileOuterRadius) { continue; }
00124
00125
00126
00127 aPrim.StartLoop();
00128 while((pNucleon = aPrim.GetNextNucleon()))
00129 {
00130 if(pNucleon->GetPosition().y()>projectileHorizon)
00131 {
00132
00133 ++particlesFromProjectile;
00134 if(pNucleon->GetParticleType() == proton)
00135 {
00136 ++chargedFromProjectile;
00137 }
00138 }
00139 }
00140 }
00141 ++hits;
00142
00143
00144 G4double targetHorizon = impactParameter-projectileOuterRadius;
00145 G4int chargedFromTarget = 0;
00146 G4int particlesFromTarget = 0;
00147 aTarg.StartLoop();
00148 while((pNucleon = aTarg.GetNextNucleon()))
00149 {
00150 if(pNucleon->GetPosition().y()>targetHorizon)
00151 {
00152
00153 ++particlesFromTarget;
00154 if(pNucleon->GetParticleType() == proton)
00155 {
00156 ++chargedFromTarget;
00157 }
00158 }
00159 }
00160
00161
00162
00163 G4ThreeVector momentum = thePrimary.Get4Momentum().vect();
00164 G4double w = (G4double)particlesFromProjectile/(G4double)aProjectileA;
00165
00166 G4double projTotEnergy = thePrimary.GetTotalEnergy();
00167 G4double targetMass = G4NucleiProperties::GetNuclearMass(aTargetA, aTargetZ);
00168 G4LorentzVector fragment4Momentum(momentum*w, projTotEnergy*w + targetMass);
00169
00170
00171 G4Fragment anInitialState(aTargetA+particlesFromProjectile,
00172 aTargetZ+chargedFromProjectile,
00173 fragment4Momentum);
00174
00175
00176 anInitialState.SetNumberOfExcitedParticle(particlesFromProjectile+particlesFromTarget,
00177 chargedFromProjectile + chargedFromTarget);
00178 anInitialState.SetNumberOfHoles(particlesFromProjectile+particlesFromTarget,
00179 chargedFromProjectile + chargedFromTarget);
00180 G4double time = thePrimary.GetGlobalTime();
00181 anInitialState.SetCreationTime(time);
00182
00183
00184 G4ReactionProductVector* thePreCompoundResult = theModel->DeExcite(anInitialState);
00185
00186
00187 G4ReactionProductVector * theExcitationResult = 0;
00188 if(particlesFromProjectile < aProjectileA)
00189 {
00190 G4LorentzVector residual4Momentum(momentum*(1.0-w), projTotEnergy*(1.0-w));
00191
00192 G4Fragment initialState2(aProjectileA-particlesFromProjectile,
00193 aProjectileZ-chargedFromProjectile,
00194 residual4Momentum );
00195
00196
00197 G4int pinit = (aProjectileA-particlesFromProjectile)/2;
00198 G4int cinit = (aProjectileZ-chargedFromProjectile)/2;
00199
00200 initialState2.SetNumberOfExcitedParticle(pinit,cinit);
00201 initialState2.SetNumberOfHoles(pinit,cinit);
00202 initialState2.SetCreationTime(time);
00203
00204 theExcitationResult = theHandler->BreakItUp(initialState2);
00205 }
00206
00207
00208 G4int nexc = 0;
00209 G4int npre = 0;
00210 if(theExcitationResult) { nexc = theExcitationResult->size(); }
00211 if(thePreCompoundResult) { npre = thePreCompoundResult->size();}
00212
00213 if(nexc > 0) {
00214 for(G4int k=0; k<nexc; ++k) {
00215 G4ReactionProduct* p = (*theExcitationResult)[k];
00216 theResult.AddSecondary(new G4DynamicParticle(p->GetDefinition(),p->GetMomentum()));
00217 delete p;
00218 }
00219 }
00220
00221 if(npre > 0) {
00222 for(G4int k=0; k<npre; ++k) {
00223 G4ReactionProduct* p = (*thePreCompoundResult)[k];
00224 theResult.AddSecondary(new G4DynamicParticle(p->GetDefinition(),p->GetMomentum()));
00225 delete p;
00226 }
00227 }
00228
00229 delete thePreCompoundResult;
00230 delete theExcitationResult;
00231
00232
00233 return &theResult;
00234
00235 }