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 #include <vector>
00048
00049 #include "G4CascadeRecoilMaker.hh"
00050 #include "globals.hh"
00051 #include "G4SystemOfUnits.hh"
00052 #include "G4CascadParticle.hh"
00053 #include "G4CascadeCheckBalance.hh"
00054 #include "G4CollisionOutput.hh"
00055 #include "G4Fragment.hh"
00056 #include "G4InuclElementaryParticle.hh"
00057 #include "G4InuclNuclei.hh"
00058 #include "G4InuclParticle.hh"
00059 #include "G4InuclSpecialFunctions.hh"
00060 #include "G4LorentzVector.hh"
00061
00062 using namespace G4InuclSpecialFunctions;
00063
00064
00065
00066
00067 G4CascadeRecoilMaker::G4CascadeRecoilMaker(G4double tolerance)
00068 : G4VCascadeCollider("G4CascadeRecoilMaker"),
00069 excTolerance(tolerance), inputEkin(0.),
00070 recoilA(0), recoilZ(0), excitationEnergy(0.) {
00071 balance = new G4CascadeCheckBalance(tolerance, tolerance, theName);
00072 }
00073
00074 G4CascadeRecoilMaker::~G4CascadeRecoilMaker() {
00075 delete balance;
00076 }
00077
00078
00079
00080
00081 void G4CascadeRecoilMaker::collide(G4InuclParticle* bullet,
00082 G4InuclParticle* target,
00083 G4CollisionOutput& output) {
00084 if (verboseLevel > 1)
00085 G4cout << " >>> G4CascadeRecoilMaker::collide" << G4endl;
00086
00087
00088 inputEkin = bullet ? bullet->getKineticEnergy() : 0.;
00089
00090 balance->setVerboseLevel(verboseLevel);
00091 balance->collide(bullet, target, output);
00092 fillRecoil();
00093 }
00094
00095
00096
00097 void G4CascadeRecoilMaker::collide(G4InuclParticle* bullet,
00098 G4InuclParticle* target,
00099 G4CollisionOutput& output,
00100 const std::vector<G4CascadParticle>& cparticles) {
00101 if (verboseLevel > 1)
00102 G4cout << " >>> G4CascadeRecoilMaker::collide(<EP>,<CP>)" << G4endl;
00103
00104
00105 inputEkin = bullet ? bullet->getKineticEnergy() : 0.;
00106
00107 balance->setVerboseLevel(verboseLevel);
00108 balance->collide(bullet, target, output, cparticles);
00109 fillRecoil();
00110 }
00111
00112
00113
00114
00115
00116 void G4CascadeRecoilMaker::fillRecoil() {
00117 recoilZ = -(balance->deltaQ());
00118 recoilA = -(balance->deltaB());
00119 recoilMomentum = -(balance->deltaLV());
00120
00121 theExcitons.clear();
00122
00123
00124 if (!goodFragment()) excitationEnergy = 0.;
00125 else excitationEnergy = deltaM() * GeV;
00126
00127
00128 if (std::abs(excitationEnergy) < excTolerance) excitationEnergy = 0.;
00129
00130 if (verboseLevel > 2) {
00131 G4cout << " recoil px " << recoilMomentum.px()
00132 << " py " << recoilMomentum.py() << " pz " << recoilMomentum.pz()
00133 << " E " << recoilMomentum.e() << " baryon " << recoilA
00134 << " charge " << recoilZ
00135 << "\n recoil mass " << recoilMomentum.m()
00136 << " 'excitation' energy " << excitationEnergy << G4endl;
00137 }
00138 }
00139
00140
00141
00142
00143 G4InuclNuclei*
00144 G4CascadeRecoilMaker::makeRecoilNuclei(G4InuclParticle::Model model) {
00145 if (verboseLevel > 1)
00146 G4cout << " >>> G4CascadeRecoilMaker::makeRecoilNuclei" << G4endl;
00147
00148 if (!goodRecoil()) {
00149 if (verboseLevel > 2 && !wholeEvent())
00150 G4cout << theName << ": event recoil is not a physical nucleus" << G4endl;
00151
00152 return 0;
00153 }
00154
00155 theRecoilNuclei.fill(recoilMomentum, recoilA, recoilZ,
00156 excitationEnergy, model);
00157 theRecoilNuclei.setExitonConfiguration(theExcitons);
00158
00159 return &theRecoilNuclei;
00160 }
00161
00162
00163
00164
00165 G4Fragment* G4CascadeRecoilMaker::makeRecoilFragment() {
00166 if (verboseLevel > 1)
00167 G4cout << " >>> G4CascadeRecoilMaker::makeRecoilFragment" << G4endl;
00168
00169 if (!goodRecoil()) {
00170 if (verboseLevel > 2 && !wholeEvent())
00171 G4cout << theName << ": event recoil is not a physical nucleus" << G4endl;
00172
00173 return 0;
00174 }
00175
00176 theRecoilFragment.SetZandA_asInt(recoilZ, recoilA);
00177
00178
00179 G4double fragMass =
00180 G4InuclNuclei::getNucleiMass(recoilA,recoilZ) + excitationEnergy/GeV;
00181
00182 G4LorentzVector fragMom; fragMom.setVectM(recoilMomentum.vect(), fragMass);
00183 theRecoilFragment.SetMomentum(fragMom*GeV);
00184
00185
00186
00187 G4int nholes = theExcitons.protonHoles+theExcitons.neutronHoles;
00188 theRecoilFragment.SetNumberOfHoles(nholes, theExcitons.protonHoles);
00189
00190 G4int nexcit = (theExcitons.protonQuasiParticles
00191 + theExcitons.neutronQuasiParticles);
00192 theRecoilFragment.SetNumberOfExcitedParticle(nexcit,
00193 theExcitons.protonQuasiParticles);
00194
00195 return &theRecoilFragment;
00196 }
00197
00198
00199
00200
00201 G4double G4CascadeRecoilMaker::deltaM() const {
00202 G4double nucMass = G4InuclNuclei::getNucleiMass(recoilA,recoilZ);
00203 return (recoilMomentum.m() - nucMass);
00204 }
00205
00206
00207
00208
00209 G4bool G4CascadeRecoilMaker::goodFragment() const {
00210 return (recoilA>0 && recoilZ>=0 && recoilA >= recoilZ);
00211 }
00212
00213 G4bool G4CascadeRecoilMaker::goodRecoil() const {
00214 return (goodFragment() && excitationEnergy > -excTolerance);
00215 }
00216
00217 G4bool G4CascadeRecoilMaker::wholeEvent() const {
00218 if (verboseLevel > 2) {
00219 G4cout << " >>> G4CascadeRecoilMaker::wholeEvent:"
00220 << " A " << recoilA << " Z " << recoilZ
00221 << " P " << recoilMomentum.rho() << " E " << recoilMomentum.e()
00222 << "\n wholeEvent returns "
00223 << (recoilA==0 && recoilZ==0 &&
00224 recoilMomentum.rho() < excTolerance/GeV &&
00225 std::abs(recoilMomentum.e()) < excTolerance/GeV) << G4endl;
00226 }
00227
00228 return (recoilA==0 && recoilZ==0 &&
00229 recoilMomentum.rho() < excTolerance/GeV &&
00230 std::abs(recoilMomentum.e()) < excTolerance/GeV);
00231 }
00232
00233
00234
00235 G4bool G4CascadeRecoilMaker::goodNucleus() const {
00236 if (verboseLevel > 2) {
00237 G4cout << " >>> G4CascadeRecoilMaker::goodNucleus" << G4endl;
00238 }
00239
00240 const G4double minExcitation = 0.1*keV;
00241 const G4double reasonableExcitation = 7.0;
00242 const G4double fractionalExcitation = 0.2;
00243
00244 if (!goodRecoil()) {
00245 if (verboseLevel>2) {
00246 if (!goodFragment()) G4cerr << " goodNucleus: invalid A/Z" << G4endl;
00247 else if (excitationEnergy < -excTolerance)
00248 G4cerr << " goodNucleus: negative excitation" << G4endl;
00249 }
00250 return false;
00251 }
00252
00253 if (excitationEnergy <= minExcitation) return true;
00254
00255
00256 G4double dm = bindingEnergy(recoilA,recoilZ);
00257 G4double exc_max0z = fractionalExcitation * inputEkin*GeV;
00258 G4double exc_dm = reasonableExcitation * dm;
00259 G4double exc_max = (exc_max0z > exc_dm) ? exc_max0z : exc_dm;
00260
00261 if (verboseLevel > 3) {
00262 G4cout << " eexs " << excitationEnergy << " max " << exc_max
00263 << " dm " << dm << G4endl;
00264 }
00265
00266 if (verboseLevel > 2 && excitationEnergy >= exc_max)
00267 G4cerr << " goodNucleus: too much excitation" << G4endl;
00268
00269 return (excitationEnergy < exc_max);
00270 }