102 theParticleChange.Clear();
103 const G4HadProjectile* aParticle = &aTrack;
104 G4double ekin = aParticle->GetKineticEnergy();
106 G4int
A = targetNucleus.GetA_asInt();
107 G4int
Z = targetNucleus.GetZ_asInt();
109 G4double m_K0S = G4KaonZeroShort::KaonZeroShort()->GetPDGMass();
110 G4double m_L = G4AntiLambda::AntiLambda()->GetPDGMass();
128 G4LorentzVector lv1 = aParticle->Get4Momentum();
129 edm::LogVerbatim(
"CMSSWNeutronAnnih") <<
"The neutron Fermi momentum (mag, x, y, z) " 130 << targetNucleus.GetFermiMomentum().mag() / MeV <<
" " 131 << targetNucleus.GetFermiMomentum().x() / MeV <<
" " 132 << targetNucleus.GetFermiMomentum().y() / MeV <<
" " 133 << targetNucleus.GetFermiMomentum().z() / MeV;
137 G4double k_neutron =
momDistr(G4UniformRand());
138 G4double momentum_neutron = 0.1973 * GeV * k_neutron;
140 G4double theta_neutron = TMath::ACos(2 * G4UniformRand() - 1);
141 G4double phi_neutron = 2. *
TMath::Pi() * G4UniformRand();
143 G4double p_neutron_x = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Cos(phi_neutron);
144 G4double p_neutron_y = momentum_neutron * TMath::Sin(theta_neutron) * TMath::Sin(phi_neutron);
145 G4double p_neutron_z = momentum_neutron * TMath::Cos(theta_neutron);
148 G4LorentzVector lv0(p_neutron_x,
151 sqrt(
pow(G4Neutron::Neutron()->GetPDGMass(), 2) + momentum_neutron * momentum_neutron));
154 G4double BENeutronInNucleus = 0;
156 BENeutronInNucleus = G4NucleiProperties::GetBindingEnergy(
A,
Z) / (
A);
158 edm::LogVerbatim(
"CMSSWNeutronAnnih") <<
"BE of nucleon in the nucleus (GeV): " << BENeutronInNucleus / GeV;
160 G4LorentzVector lvBE(0, 0, 0, BENeutronInNucleus / GeV);
161 G4LorentzVector lv = lv0 + lv1 - lvBE;
164 G4double etot = lv0.e() + lv1.e() - lvBE.e();
165 if (etot < theK0S->GetPDGMass() +
theAntiL->GetPDGMass()) {
166 theParticleChange.SetEnergyChange(ekin);
167 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
168 return &theParticleChange;
171 float newIonMass = targetNucleus.AtomicMass(
A - 1,
Z) * 931.5 * MeV;
173 G4LorentzVector nlvIon(0, 0, 0, newIonMass);
175 G4double theta_KS0_star = TMath::ACos(2 * G4UniformRand() - 1);
176 G4double phi_KS0_star = 2. *
TMath::Pi() * G4UniformRand();
178 G4double p_K0S_star_x = TMath::Sin(theta_KS0_star) * TMath::Cos(phi_KS0_star);
179 G4double p_K0S_star_y = TMath::Sin(theta_KS0_star) * TMath::Sin(phi_KS0_star);
180 G4double p_K0S_star_z = TMath::Cos(theta_KS0_star);
182 G4ThreeVector
p(p_K0S_star_x, p_K0S_star_y, p_K0S_star_z);
184 double m0_2 =
m0 *
m0;
185 double m1_2 = m_K0S * m_K0S;
186 double m2_2 = m_L * m_L;
188 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);
189 double p2 =
p.mag2();
191 G4LorentzVector nlvK0S(
p,
sqrt(
p2 + m1_2));
192 G4LorentzVector nlvAntiL(-
p,
sqrt(
p2 + m2_2));
195 nlvK0S.boost(lv.boostVector());
196 nlvAntiL.boost(lv.boostVector());
199 theParticleChange.SetStatusChange(stopAndKill);
202 G4DynamicParticle* aSec1 =
new G4DynamicParticle(
theK0S, nlvK0S);
203 theParticleChange.AddSecondary(aSec1);
204 G4DynamicParticle* aSec2 =
new G4DynamicParticle(
theAntiL, nlvAntiL);
205 theParticleChange.AddSecondary(aSec2);
207 const G4ParticleDefinition* theRemainingNucleusDef =
theProton;
209 theRemainingNucleusDef = G4IonTable::GetIonTable()->GetIon(
Z,
A - 1);
210 G4DynamicParticle* aSec3 =
new G4DynamicParticle(theRemainingNucleusDef, nlvIon);
211 theParticleChange.AddSecondary(aSec3);
214 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)