105 theParticleChange.Clear();
106 const G4HadProjectile* aParticle = &aTrack;
107 G4double ekin = aParticle->GetKineticEnergy();
109 G4int
A = targetNucleus.GetA_asInt();
110 G4int
Z = targetNucleus.GetZ_asInt();
112 G4double m_K0S = G4KaonZeroShort::KaonZeroShort()->GetPDGMass();
113 G4double m_L = G4AntiLambda::AntiLambda()->GetPDGMass();
131 G4LorentzVector lv1 = aParticle->Get4Momentum();
132 edm::LogVerbatim(
"CMSSWNeutronAnnih") <<
"The neutron Fermi momentum (mag, x, y, z) " 133 << targetNucleus.GetFermiMomentum().mag() / MeV <<
" " 134 << targetNucleus.GetFermiMomentum().x() / MeV <<
" " 135 << targetNucleus.GetFermiMomentum().y() / MeV <<
" " 136 << targetNucleus.GetFermiMomentum().z() / MeV;
140 G4double k_neutron =
momDistr(G4UniformRand());
141 G4double momentum_neutron = 0.1973 * GeV * k_neutron;
143 G4double theta_neutron = TMath::ACos(2 * G4UniformRand() - 1);
144 G4double phi_neutron = 2. *
TMath::Pi() * G4UniformRand();
146 G4double p_neutron_x = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Cos(phi_neutron);
147 G4double p_neutron_y = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Sin(phi_neutron);
148 G4double p_neutron_z = momentum_neutron * TMath::Cos(theta_neutron);
151 G4LorentzVector lv0(p_neutron_x,
154 sqrt(
pow(G4Neutron::Neutron()->GetPDGMass(), 2) + momentum_neutron * momentum_neutron));
157 G4double BENeutronInNucleus = 0;
159 BENeutronInNucleus = G4NucleiProperties::GetBindingEnergy(
A,
Z) / (
A);
161 edm::LogVerbatim(
"CMSSWNeutronAnnih") <<
"BE of nucleon in the nucleus (GeV): " << BENeutronInNucleus / GeV;
163 G4LorentzVector lvBE(0, 0, 0, BENeutronInNucleus / GeV);
164 G4LorentzVector lv = lv0 + lv1 - lvBE;
167 G4double etot = lv0.e() + lv1.e() - lvBE.e();
168 if (etot < theK0S->GetPDGMass() +
theAntiL->GetPDGMass()) {
169 theParticleChange.SetEnergyChange(ekin);
170 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
171 return &theParticleChange;
174 float newIonMass = targetNucleus.AtomicMass(
A - 1,
Z) * 931.5 * MeV;
176 G4LorentzVector nlvIon(0, 0, 0, newIonMass);
178 G4double theta_KS0_star = TMath::ACos(2 * G4UniformRand() - 1);
179 G4double phi_KS0_star = 2. *
TMath::Pi() * G4UniformRand();
181 G4double p_K0S_star_x = TMath::Sin(theta_KS0_star) * TMath::Cos(phi_KS0_star);
182 G4double p_K0S_star_y = TMath::Sin(theta_KS0_star) * TMath::Sin(phi_KS0_star);
183 G4double p_K0S_star_z = TMath::Cos(theta_KS0_star);
185 G4ThreeVector
p(p_K0S_star_x, p_K0S_star_y, p_K0S_star_z);
187 double m0_2 =
m0 *
m0;
188 double m1_2 = m_K0S * m_K0S;
189 double m2_2 = m_L * m_L;
191 p *= 0.5 /
m0 *
sqrt(m0_2 * m0_2 + m1_2 * m1_2 + m2_2 * m2_2 - 2 * m0_2 * m1_2 - 2 * m0_2 * m2_2 - 2 * m1_2 * m2_2);
192 double p2 =
p.mag2();
194 G4LorentzVector nlvK0S(
p,
sqrt(
p2 + m1_2));
195 G4LorentzVector nlvAntiL(-
p,
sqrt(
p2 + m2_2));
198 nlvK0S.boost(lv.boostVector());
199 nlvAntiL.boost(lv.boostVector());
202 theParticleChange.SetStatusChange(stopAndKill);
205 G4DynamicParticle* aSec1 =
new G4DynamicParticle(
theK0S, nlvK0S);
206 theParticleChange.AddSecondary(aSec1);
207 G4DynamicParticle* aSec2 =
new G4DynamicParticle(
theAntiL, nlvAntiL);
208 theParticleChange.AddSecondary(aSec2);
210 const G4ParticleDefinition* theRemainingNucleusDef =
theProton;
212 theRemainingNucleusDef = G4IonTable::GetIonTable()->GetIon(
Z,
A - 1);
213 G4DynamicParticle* aSec3 =
new G4DynamicParticle(theRemainingNucleusDef, nlvIon);
214 theParticleChange.AddSecondary(aSec3);
217 return &theParticleChange;
Log< level::Info, true > LogVerbatim
G4ParticleDefinition * theProton
G4double momDistr(G4double x_in)
G4ParticleDefinition * theAntiL
G4ParticleDefinition * theK0S
Power< A, B >::type pow(const A &a, const B &b)