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 #include "G4PiMinusStopAbsorption.hh"
00036 #include <vector>
00037
00038 #include "globals.hh"
00039 #include "Randomize.hh"
00040 #include "G4NucleiProperties.hh"
00041 #include "G4ParticleTypes.hh"
00042 #include "G4Nucleus.hh"
00043 #include "G4ReactionKinematics.hh"
00044 #include "G4DynamicParticleVector.hh"
00045 #include "G4ParticleDefinition.hh"
00046 #include "G4Proton.hh"
00047 #include "G4Neutron.hh"
00048 #include "G4ThreeVector.hh"
00049 #include "G4HadronicDeprecate.hh"
00050
00051
00052
00053 G4PiMinusStopAbsorption::G4PiMinusStopAbsorption(G4PiMinusStopMaterial* materialAlgo,
00054 const G4double Z, const G4double A)
00055
00056 {
00057 G4HadronicDeprecate("G4PiMinusStopAbsorption");
00058 _materialAlgo = materialAlgo;
00059 _nucleusZ = Z;
00060 _nucleusA = A;
00061 _level = 0;
00062 _absorptionProducts = new G4DynamicParticleVector();
00063 }
00064
00065
00066
00067 G4PiMinusStopAbsorption::~G4PiMinusStopAbsorption()
00068 {
00069
00070 delete _materialAlgo;
00071
00072
00073 delete _absorptionProducts;
00074 }
00075
00076 G4DynamicParticleVector* G4PiMinusStopAbsorption::DoAbsorption()
00077 {
00078 std::vector<G4ParticleDefinition*>* defNucleons = _materialAlgo->DefinitionVector();
00079
00080 G4double newA = _nucleusA;
00081 G4double newZ = _nucleusZ;
00082
00083 if (defNucleons != 0)
00084 {
00085 for (unsigned int i=0; i<defNucleons->size(); i++)
00086 {
00087 if ( (*defNucleons)[i] == G4Proton::Proton())
00088 {
00089 newA = newA - 1;
00090 newZ = newZ - 1;
00091 }
00092 if ((*defNucleons)[i] == G4Neutron::Neutron())
00093 { newA = newA - 1; }
00094 }
00095 }
00096
00097 G4double binding = G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(_nucleusA) ,static_cast<G4int>(_nucleusZ)) / _nucleusA;
00098 G4double mass = G4NucleiProperties::GetNuclearMass(static_cast<G4int>(newA),static_cast<G4int>(newZ));
00099
00100
00101 std::vector<G4LorentzVector*>* p4Nucleons = _materialAlgo->P4Vector(binding,mass);
00102
00103 if (defNucleons != 0 && p4Nucleons != 0)
00104 {
00105 unsigned int nNucleons = p4Nucleons->size();
00106
00107 G4double seen = _materialAlgo->FinalNucleons() / 2.;
00108 G4int maxN = nNucleons;
00109 if (defNucleons->size() < nNucleons) { maxN = defNucleons->size(); }
00110
00111 for (G4int i=0; i<maxN; i++)
00112 {
00113 G4DynamicParticle* product;
00114 if ((*defNucleons)[i] == G4Proton::Proton())
00115 { product = new G4DynamicParticle(G4Proton::Proton(),*((*p4Nucleons)[i])); }
00116 else
00117 { product = new G4DynamicParticle(G4Neutron::Neutron(),*((*p4Nucleons)[i])); }
00118 G4double ranflat = G4UniformRand();
00119
00120 if (ranflat < seen)
00121 { _absorptionProducts->push_back(product); }
00122 else
00123 { delete product; }
00124 }
00125 }
00126
00127 return _absorptionProducts;
00128
00129 }
00130
00131 G4ThreeVector G4PiMinusStopAbsorption::RecoilMomentum()
00132 {
00133 G4ThreeVector pProducts(0.,0.,0.);
00134
00135 for (unsigned int i = 0; i< _absorptionProducts->size(); i++)
00136 {
00137 pProducts = pProducts + (*_absorptionProducts)[i]->GetMomentum();
00138 }
00139 return pProducts;
00140 }
00141
00142
00143 G4int G4PiMinusStopAbsorption::NProtons()
00144 {
00145 G4int n = 0;
00146 G4int entries = _absorptionProducts->size();
00147 for (int i = 0; i<entries; i++)
00148 {
00149 if ((*_absorptionProducts)[i]->GetDefinition() == G4Proton::Proton())
00150 { n = n + 1; }
00151 }
00152 return n;
00153 }
00154
00155
00156 G4int G4PiMinusStopAbsorption::NNeutrons()
00157 {
00158 G4int n = 0;
00159 G4int entries = _absorptionProducts->size();
00160 for (int i = 0; i<entries; i++)
00161 {
00162 if ((*_absorptionProducts)[i]->GetDefinition() == G4Neutron::Neutron())
00163 { n = n + 1; }
00164 }
00165 return n;
00166 }
00167
00168
00169 G4double G4PiMinusStopAbsorption::Energy()
00170 {
00171 G4double energy = 0.;
00172 G4double productEnergy = 0.;
00173 G4ThreeVector pProducts(0.,0.,0.);
00174 G4int nN = 0;
00175 G4int nP = 0;
00176
00177
00178 G4int nAbsorptionProducts = _absorptionProducts->size();
00179
00180 for (int i = 0; i<nAbsorptionProducts; i++)
00181 {
00182 productEnergy += (*_absorptionProducts)[i]->GetKineticEnergy();
00183 pProducts = pProducts + (*_absorptionProducts)[i]->GetMomentum();
00184 if ((*_absorptionProducts)[i]->GetDefinition() == G4Neutron::Neutron()) nN++;
00185 if ((*_absorptionProducts)[i]->GetDefinition() == G4Proton::Proton()) nP++;
00186 }
00187
00188 G4double productBinding = (G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(_nucleusA),static_cast<G4int>(_nucleusZ)) / _nucleusA) * nAbsorptionProducts;
00189 G4double mass = G4NucleiProperties::GetNuclearMass(_nucleusA - (nP + nN),_nucleusZ - nP);
00190 G4double pNucleus = pProducts.mag();
00191 G4double eNucleus = std::sqrt(pNucleus*pNucleus + mass*mass);
00192 G4double tNucleus = eNucleus - mass;
00193 G4double temp =
00194 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(_nucleusA - (nP + nN)),static_cast<G4int>(_nucleusZ - nP)) -
00195 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(_nucleusA),static_cast<G4int>(_nucleusZ));
00196 energy = productEnergy + productBinding + tNucleus;
00197
00198 if (_level > 0)
00199 {
00200 std::cout << "E products " << productEnergy
00201 << " Binding " << productBinding << " " << temp << " "
00202 << " Tnucleus " << tNucleus
00203 << " energy = " << energy << G4endl;
00204 }
00205
00206 return energy;
00207 }
00208
00209 void G4PiMinusStopAbsorption::SetVerboseLevel(G4int level)
00210 {
00211 _level = level;
00212 return;
00213 }