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 #include <vector>
00035
00036 #include "G4PiMinusStopMaterial.hh"
00037
00038 #include "globals.hh"
00039 #include "G4ios.hh"
00040 #include "Randomize.hh"
00041 #include "G4PhysicalConstants.hh"
00042 #include "G4Proton.hh"
00043 #include "G4Neutron.hh"
00044 #include "G4PionMinus.hh"
00045 #include "G4ParticleTypes.hh"
00046 #include "G4ReactionKinematics.hh"
00047 #include "G4DynamicParticleVector.hh"
00048 #include "G4LorentzVector.hh"
00049 #include "G4PiMinusStopMaterial.hh"
00050 #include "G4DistributionGenerator.hh"
00051
00052
00053
00054
00055 G4PiMinusStopMaterial::G4PiMinusStopMaterial()
00056 {
00057 _definitions = 0;
00058 _momenta = 0;
00059 _distributionE = 0;
00060 _distributionAngle = 0;
00061 theR = 0.5;
00062 }
00063
00064
00065
00066
00067 G4PiMinusStopMaterial::~G4PiMinusStopMaterial()
00068 {
00069 if (_definitions != 0) delete _definitions;
00070 _definitions = 0;
00071
00072
00073 if (_momenta != 0) {
00074 for (unsigned int i=0; i<_momenta->size(); i++) delete(*_momenta)[i];
00075 delete _momenta;
00076 }
00077
00078 delete _distributionE;
00079 delete _distributionAngle;
00080 }
00081
00082 std::vector<G4ParticleDefinition*>* G4PiMinusStopMaterial::DefinitionVector()
00083 {
00084
00085 _definitions->push_back(G4Neutron::Neutron());
00086
00087 G4double ranflat = G4UniformRand();
00088 if (ranflat < theR)
00089 { _definitions->push_back(G4Proton::Proton()); }
00090 else
00091 { _definitions->push_back(G4Neutron::Neutron()); }
00092
00093 return _definitions;
00094
00095 }
00096
00097 std::vector<G4LorentzVector*>*
00098 G4PiMinusStopMaterial::P4Vector(const G4double binding,
00099 const G4double massNucleus)
00100 {
00101
00102
00103
00104
00105 G4double eKin1;
00106 G4double eKin2;
00107 G4double eRecoil;
00108
00109
00110 G4int nNucleons = 2;
00111 G4double availableE = G4PionMinus::PionMinus()->GetPDGMass() - nNucleons * binding;
00112 G4LorentzVector p1;
00113 G4LorentzVector p2;
00114
00115 do
00116 {
00117 G4double ranflat;
00118 G4double p;
00119 G4double energy;
00120 G4double mass;
00121
00122 ranflat = G4UniformRand();
00123 eKin1 = _distributionE->Generate(ranflat);
00124 mass = (*_definitions)[0]->GetPDGMass();
00125 energy = eKin1 + mass;
00126 p = std::sqrt(energy*energy - mass*mass);
00127 G4double theta1 = pi*G4UniformRand();
00128 G4double phi1 = GenerateAngle(2.*pi);
00129 p1 = MakeP4(p,theta1,phi1,energy);
00130
00131 ranflat = G4UniformRand();
00132 eKin2 = _distributionE->Generate(ranflat);
00133 mass = (*_definitions)[1]->GetPDGMass();
00134 energy = eKin2 + mass;
00135 p = std::sqrt(energy*energy - mass*mass);
00136 ranflat = G4UniformRand();
00137 G4double opAngle = _distributionAngle->Generate(ranflat);
00138 G4double theta2 = theta1 + opAngle;
00139 G4double phi2 = phi1 + opAngle;
00140
00141 p2 = MakeP4(p,theta2,phi2,energy);
00142
00143 G4double pNucleus = (p1.vect() + p2.vect()).mag();
00144 eRecoil = std::sqrt(pNucleus*pNucleus + massNucleus*massNucleus) - massNucleus;
00145
00146
00147
00148
00149
00150
00151
00152
00153 } while ((eKin1 + eKin2 + eRecoil) > availableE);
00154
00155
00156 if (_momenta != 0) {
00157 _momenta->push_back(new G4LorentzVector(p1));
00158 _momenta->push_back(new G4LorentzVector(p2));
00159 }
00160
00161 return _momenta;
00162
00163 }
00164
00165 G4double G4PiMinusStopMaterial::GenerateAngle(G4double x)
00166 {
00167 G4double ranflat = G4UniformRand();
00168 G4double value = ranflat * x;
00169 return value;
00170 }
00171
00172 G4LorentzVector G4PiMinusStopMaterial::MakeP4(G4double p, G4double theta, G4double phi, G4double e)
00173 {
00174
00175 G4double px = p * std::sin(theta) * std::cos(phi);
00176 G4double py = p * std::sin(theta) * std::sin(phi);
00177 G4double pz = p * std::cos(theta);
00178 G4LorentzVector p4(px,py,pz,e);
00179 return p4;
00180 }
00181
00182 G4double G4PiMinusStopMaterial::RecoilEnergy(const G4double mass)
00183 {
00184 G4ThreeVector p(0.,0.,0.);
00185
00186
00187 if (_momenta != 0) {
00188 for (unsigned int i = 0; i< _momenta->size(); i++)
00189 {
00190 p = p + (*_momenta)[i]->vect();
00191 }
00192 }
00193
00194 G4double pNucleus = p.mag();
00195 G4double eNucleus = std::sqrt(pNucleus*pNucleus + mass*mass);
00196
00197 return eNucleus;
00198 }
00199