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
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 #ifndef G4NUCLEI_MODEL_HH
00068 #define G4NUCLEI_MODEL_HH
00069
00070 #include <algorithm>
00071 #include <vector>
00072 #include <CLHEP/Units/SystemOfUnits.h>
00073
00074 #include "G4InuclElementaryParticle.hh"
00075 #include "G4CascadParticle.hh"
00076 #include "G4CollisionOutput.hh"
00077 #include "G4LorentzConvertor.hh"
00078
00079 class G4InuclNuclei;
00080 class G4ElementaryParticleCollider;
00081
00082 class G4NucleiModel {
00083 public:
00084 G4NucleiModel();
00085 G4NucleiModel(G4int a, G4int z);
00086 explicit G4NucleiModel(G4InuclNuclei* nuclei);
00087
00088 ~G4NucleiModel();
00089
00090 void setVerboseLevel(G4int verbose) { verboseLevel = verbose; }
00091
00092 void generateModel(G4InuclNuclei* nuclei);
00093 void generateModel(G4int a, G4int z);
00094
00095
00096 void reset(G4int nHitNeutrons=0, G4int nHitProtons=0,
00097 const std::vector<G4ThreeVector>* hitPoints=0);
00098
00099 void printModel() const;
00100
00101 G4double getDensity(G4int ip, G4int izone) const {
00102 return nucleon_densities[ip - 1][izone];
00103 }
00104
00105 G4double getFermiMomentum(G4int ip, G4int izone) const {
00106 return fermi_momenta[ip - 1][izone];
00107 }
00108
00109 G4double getFermiKinetic(G4int ip, G4int izone) const;
00110
00111 G4double getPotential(G4int ip, G4int izone) const {
00112 if (ip == 9) return 0.0;
00113 G4int ip0 = ip < 3 ? ip - 1 : 2;
00114 if (ip > 10 && ip < 18) ip0 = 3;
00115 if (ip > 20) ip0 = 4;
00116 return izone < number_of_zones ? zone_potentials[ip0][izone] : 0.0;
00117 }
00118
00119
00120 G4double getRadiusUnits() const { return radiusUnits*CLHEP::fermi; }
00121
00122 G4double getRadius() const { return nuclei_radius; }
00123 G4double getRadius(G4int izone) const {
00124 return ( (izone<0) ? 0.
00125 : (izone<number_of_zones) ? zone_radii[izone] : nuclei_radius);
00126 }
00127 G4double getVolume(G4int izone) const {
00128 return ( (izone<0) ? 0.
00129 : (izone<number_of_zones) ? zone_volumes[izone] : nuclei_volume);
00130 }
00131
00132 G4int getNumberOfZones() const { return number_of_zones; }
00133 G4int getZone(G4double r) const {
00134 for (G4int iz=0; iz<number_of_zones; iz++) if (r<zone_radii[iz]) return iz;
00135 return number_of_zones;
00136 }
00137
00138 G4int getNumberOfNeutrons() const { return neutronNumberCurrent; }
00139 G4int getNumberOfProtons() const { return protonNumberCurrent; }
00140
00141 G4bool empty() const {
00142 return neutronNumberCurrent < 1 && protonNumberCurrent < 1;
00143 }
00144
00145 G4bool stillInside(const G4CascadParticle& cparticle) {
00146 return cparticle.getCurrentZone() < number_of_zones;
00147 }
00148
00149
00150 G4CascadParticle initializeCascad(G4InuclElementaryParticle* particle);
00151
00152 typedef std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > modelLists;
00153
00154 void initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target,
00155 modelLists& output);
00156
00157
00158 std::pair<G4int, G4int> getTypesOfNucleonsInvolved() const {
00159 return std::pair<G4int, G4int>(current_nucl1, current_nucl2);
00160 }
00161
00162 void generateParticleFate(G4CascadParticle& cparticle,
00163 G4ElementaryParticleCollider* theEPCollider,
00164 std::vector<G4CascadParticle>& cascade);
00165
00166 G4bool isProjectile(const G4CascadParticle& cparticle) const;
00167 G4bool worthToPropagate(const G4CascadParticle& cparticle) const;
00168
00169 G4InuclElementaryParticle generateNucleon(G4int type, G4int zone) const;
00170
00171 G4LorentzVector generateNucleonMomentum(G4int type, G4int zone) const;
00172
00173 G4double absorptionCrossSection(G4double e, G4int type) const;
00174 G4double totalCrossSection(G4double ke, G4int rtype) const;
00175
00176 private:
00177 G4int verboseLevel;
00178
00179 G4bool passFermi(const std::vector<G4InuclElementaryParticle>& particles,
00180 G4int zone);
00181
00182 G4bool passTrailing(const G4ThreeVector& hit_position);
00183
00184 void boundaryTransition(G4CascadParticle& cparticle);
00185
00186 G4InuclElementaryParticle generateQuasiDeutron(G4int type1,
00187 G4int type2,
00188 G4int zone) const;
00189
00190 typedef std::pair<G4InuclElementaryParticle, G4double> partner;
00191
00192 std::vector<partner> thePartners;
00193 void generateInteractionPartners(G4CascadParticle& cparticle);
00194
00195
00196 static G4bool sortPartners(const partner& p1, const partner& p2) {
00197 return (p2.second > p1.second);
00198 }
00199
00200
00201 void fillBindingEnergies();
00202
00203 void fillZoneRadii(G4double nuclearRadius);
00204
00205 G4double fillZoneVolumes(G4double nuclearRadius);
00206
00207 void fillPotentials(G4int type, G4double tot_vol);
00208
00209 G4double zoneIntegralWoodsSaxon(G4double ur1, G4double ur2,
00210 G4double nuclearRadius) const;
00211
00212 G4double zoneIntegralGaussian(G4double ur1, G4double ur2,
00213 G4double nuclearRadius) const;
00214
00215 G4double getRatio(G4int ip) const;
00216
00217
00218 G4double getCurrentDensity(G4int ip, G4int izone) const;
00219
00220
00221
00222 G4double inverseMeanFreePath(const G4CascadParticle& cparticle,
00223 const G4InuclElementaryParticle& target);
00224
00225
00226
00227 G4double generateInteractionLength(const G4CascadParticle& cparticle,
00228 G4double path, G4double invmfp) const;
00229
00230
00231 G4LorentzConvertor dummy_convertor;
00232 G4CollisionOutput EPCoutput;
00233
00234 std::vector<G4InuclElementaryParticle> qdeutrons;
00235 std::vector<G4double> acsecs;
00236
00237 std::vector<G4ThreeVector> coordinates;
00238 std::vector<G4LorentzVector> momentums;
00239 std::vector<G4InuclElementaryParticle> raw_particles;
00240
00241 std::vector<G4ThreeVector> collisionPts;
00242
00243
00244 G4double ur[7];
00245 G4double v[6];
00246 G4double v1[6];
00247 std::vector<G4double> rod;
00248 std::vector<G4double> pf;
00249 std::vector<G4double> vz;
00250
00251
00252 std::vector<std::vector<G4double> > nucleon_densities;
00253 std::vector<std::vector<G4double> > zone_potentials;
00254 std::vector<std::vector<G4double> > fermi_momenta;
00255 std::vector<G4double> zone_radii;
00256 std::vector<G4double> zone_volumes;
00257 std::vector<G4double> binding_energies;
00258 G4double nuclei_radius;
00259 G4double nuclei_volume;
00260 G4int number_of_zones;
00261
00262 G4int A;
00263 G4int Z;
00264 G4InuclNuclei* theNucleus;
00265
00266 G4int neutronNumber;
00267 G4int protonNumber;
00268
00269 G4int neutronNumberCurrent;
00270 G4int protonNumberCurrent;
00271
00272 G4int current_nucl1;
00273 G4int current_nucl2;
00274
00275
00276 enum PotentialType { WoodsSaxon=0, Gaussian=1 };
00277
00278
00279 static const G4double skinDepth;
00280 static const G4double radiusScale;
00281 static const G4double radiusScale2;
00282 static const G4double radiusForSmall;
00283 static const G4double radScaleAlpha;
00284 static const G4double fermiMomentum;
00285 static const G4double R_nucleon;
00286 static const G4double alfa3[3], alfa6[6];
00287 static const G4double pion_vp;
00288 static const G4double pion_vp_small;
00289 static const G4double kaon_vp;
00290 static const G4double hyperon_vp;
00291
00292
00293 static const G4double piTimes4thirds;
00294 static const G4double crossSectionUnits;
00295 static const G4double radiusUnits;
00296 };
00297
00298 #endif // G4NUCLEI_MODEL_HH