| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4 Model: Charge and strangness exchange based on G4LightMedia model
27 // 28 May 2006 V.Ivanchenko
28 //
29 // Modified:
30 // 07-Jun-06 V.Ivanchenko fix problem of rotation of final state
31 // 25-Jul-06 V.Ivanchenko add 19 MeV low energy, below which S-wave is sampled
32 // 12-Jun-12 A.Ribon fix warnings of shadowed variables
33 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log and G4Pow
34 //
35
36 #include "G4ChargeExchange.hh"
37 #include "G4ChargeExchangeXS.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4ParticleTable.hh"
41 #include "G4ParticleDefinition.hh"
42 #include "G4IonTable.hh"
43 #include "Randomize.hh"
44 #include "G4NucleiProperties.hh"
45 #include "G4DecayTable.hh"
46 #include "G4VDecayChannel.hh"
47 #include "G4DecayProducts.hh"
48 #include "G4NistManager.hh"
49 #include "G4Fragment.hh"
50 #include "G4ExcitationHandler.hh"
51 #include "G4ReactionProductVector.hh"
52
53 #include "G4Exp.hh"
54 #include "G4Log.hh"
55 #include "G4Pow.hh"
56
57 #include "G4HadronicParameters.hh"
58 #include "G4PhysicsModelCatalog.hh"
59
60 namespace
61 {
62 constexpr G4int maxN = 1000;
63 constexpr G4double emin = 2*136.9*CLHEP::MeV;
64 }
65
66 G4ChargeExchange::G4ChargeExchange(G4ChargeExchangeXS* ptr)
67 : G4HadronicInteraction("ChargeExchange"),
68 fXSection(ptr), fXSWeightFactor(1.0)
69 {
70 lowEnergyLimit = 1.*CLHEP::MeV;
71 secID = G4PhysicsModelCatalog::GetModelID( "model_ChargeExchange" );
72 nist = G4NistManager::Instance();
73 fHandler = new G4ExcitationHandler();
74 if (nullptr != fXSection) {
75 fXSWeightFactor = 1.0/fXSection->GetCrossSectionFactor();
76 }
77 }
78
79 G4ChargeExchange::~G4ChargeExchange()
80 {
81 delete fHandler;
82 }
83
84 G4HadFinalState* G4ChargeExchange::ApplyYourself(
85 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
86 {
87 theParticleChange.Clear();
88 auto part = aTrack.GetDefinition();
89 G4double ekin = aTrack.GetKineticEnergy();
90
91 G4int A = targetNucleus.GetA_asInt();
92 G4int Z = targetNucleus.GetZ_asInt();
93
94 if (ekin <= lowEnergyLimit) {
95 return &theParticleChange;
96 }
97 theParticleChange.SetWeightChange(fXSWeightFactor);
98
99 G4int projPDG = part->GetPDGEncoding();
100
101 // for hydrogen targets and positive projectile change exchange
102 // is not possible on proton, only on deuteron
103 if (1 == Z && (211 == projPDG || 321 == projPDG)) { A = 2; }
104
105 if (verboseLevel > 1)
106 G4cout << "G4ChargeExchange for " << part->GetParticleName()
107 << " PDGcode= " << projPDG << " on nucleus Z= " << Z
108 << " A= " << A << " N= " << A - Z
109 << G4endl;
110
111 G4double mass1 = G4NucleiProperties::GetNuclearMass(A, Z);
112 G4LorentzVector lv0 = aTrack.Get4Momentum();
113 G4double etot = mass1 + lv0.e();
114
115 // select final state
116 const G4ParticleDefinition* theSecondary =
117 fXSection->SampleSecondaryType(part, Z, A);
118 G4int pdg = theSecondary->GetPDGEncoding();
119
120 // omega(782) and f2(1270)
121 G4bool isShortLived = (pdg == 223 || pdg == 225);
122
123 // atomic number of the recoil nucleus
124 if (projPDG == -211) { --Z; }
125 else if (projPDG == 211) { ++Z; }
126 else if (projPDG == -321) { --Z; }
127 else if (projPDG == 321) { ++Z; }
128 else if (projPDG == 130) {
129 if (theSecondary->GetPDGCharge() > 0.0) { --Z; }
130 else { ++Z; }
131 } else {
132 // not ready for other projectile
133 return &theParticleChange;
134 }
135
136 // recoil nucleus
137 const G4ParticleDefinition* theRecoil = nullptr;
138 if (Z == 0 && A == 1) { theRecoil = G4Neutron::Neutron(); }
139 else if (Z == 1 && A == 1) { theRecoil = G4Proton::Proton(); }
140 else if (Z == 1 && A == 2) { theRecoil = G4Deuteron::Deuteron(); }
141 else if (Z == 1 && A == 3) { theRecoil = G4Triton::Triton(); }
142 else if (Z == 2 && A == 3) { theRecoil = G4He3::He3(); }
143 else if (Z == 2 && A == 4) { theRecoil = G4Alpha::Alpha(); }
144 else if (nist->GetIsotopeAbundance(Z, A) > 0.0) {
145 theRecoil = G4ParticleTable::GetParticleTable()
146 ->GetIonTable()->GetIon(Z, A, 0.0);
147 }
148
149 // check if there is enough energy for the final state
150 // and sample mass of produced state
151 const G4double mass0 = theSecondary->GetPDGMass();
152 G4double mass3 = (nullptr == theRecoil) ?
153 G4NucleiProperties::GetNuclearMass(A, Z) : theRecoil->GetPDGMass();
154 G4double mass2 = mass0;
155 if (isShortLived &&
156 !SampleMass(mass2, theSecondary->GetPDGWidth(), etot - mass3)) {
157 return &theParticleChange;
158 }
159
160 // not possible kinematically
161 if (etot <= mass2 + mass3) {
162 return &theParticleChange;
163 }
164
165 // sample kinematics
166 G4LorentzVector lv1(0.0, 0.0, 0.0, mass1);
167 G4LorentzVector lv = lv0 + lv1;
168 G4ThreeVector bst = lv.boostVector();
169 G4double ss = lv.mag2();
170
171 // tmax = 4*momCMS^2
172 G4double e2 = ss + mass2*mass2 - mass3*mass3;
173 G4double tmax = e2*e2/ss - 4*mass2*mass2;
174
175 G4double t = SampleT(theSecondary, A, tmax);
176
177 G4double phi = G4UniformRand()*CLHEP::twopi;
178 G4double cost = 1. - 2.0*t/tmax;
179
180 if (cost > 1.0) { cost = 1.0; }
181 else if(cost < -1.0) { cost = -1.0; }
182
183 G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
184
185 if (verboseLevel>1) {
186 G4cout << " t= " << t << " tmax(GeV^2)= " << tmax/(GeV*GeV)
187 << " cos(t)=" << cost << " sin(t)=" << sint << G4endl;
188 }
189 G4double momentumCMS = 0.5*std::sqrt(tmax);
190 G4LorentzVector lv2(momentumCMS*sint*std::cos(phi),
191 momentumCMS*sint*std::sin(phi),
192 momentumCMS*cost,
193 std::sqrt(momentumCMS*momentumCMS + mass2*mass2));
194
195 // kinematics in the final state, may be a warning should be added if
196 lv2.boost(bst);
197 if (lv2.e() < mass2) {
198 lv2.setE(mass2);
199 }
200 lv -= lv2;
201 if (lv.e() < mass3) {
202 lv.setE(mass3);
203 }
204
205 // prepare secondary particles
206 theParticleChange.SetStatusChange(stopAndKill);
207 theParticleChange.SetEnergyChange(0.0);
208
209 if (!isShortLived) {
210 auto aSec = new G4DynamicParticle(theSecondary, lv2);
211 theParticleChange.AddSecondary(aSec, secID);
212 } else {
213 auto channel = theSecondary->GetDecayTable()->SelectADecayChannel(mass2);
214 auto products = channel->DecayIt(mass2);
215 G4ThreeVector bst1 = lv2.boostVector();
216 G4int N = products->entries();
217 for (G4int i=0; i<N; ++i) {
218 auto p = (*products)[i];
219 auto lvp = p->Get4Momentum();
220 lvp.boost(bst1);
221 p->Set4Momentum(lvp);
222 theParticleChange.AddSecondary(p, secID);
223 }
224 delete products;
225 }
226
227 // recoil is a stable isotope
228 if (nullptr != theRecoil) {
229 auto aRec = new G4DynamicParticle(theRecoil, lv);
230 theParticleChange.AddSecondary(aRec, secID);
231 } else {
232 // recoil is an unstable fragment
233 G4Fragment frag(A, Z, lv);
234 auto products = fHandler->BreakItUp(frag);
235 for (auto & prod : *products) {
236 auto dp = new G4DynamicParticle(prod->GetDefinition(), prod->GetMomentum());
237 theParticleChange.AddSecondary(dp, secID);
238 delete prod;
239 }
240 delete products;
241 }
242 return &theParticleChange;
243 }
244
245 G4double G4ChargeExchange::SampleT(const G4ParticleDefinition*,
246 const G4int A, const G4double tmax) const
247 {
248 G4double aa, bb, cc, dd;
249 G4Pow* g4pow = G4Pow::GetInstance();
250 if (A <= 62.) {
251 aa = g4pow->powZ(A, 1.63);
252 bb = 14.5*g4pow->powZ(A, 0.66);
253 cc = 1.4*g4pow->powZ(A, 0.33);
254 dd = 10.;
255 } else {
256 aa = g4pow->powZ(A, 1.33);
257 bb = 60.*g4pow->powZ(A, 0.33);
258 cc = 0.4*g4pow->powZ(A, 0.40);
259 dd = 10.;
260 }
261 G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb;
262 G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd;
263
264 G4double t;
265 G4double y = bb;
266 if(G4UniformRand()*(x1 + x2) < x2) y = dd;
267
268 for (G4int i=0; i<maxN; ++i) {
269 t = -G4Log(G4UniformRand())/y;
270 if (t <= tmax) { return t; }
271 }
272 return 0.0;
273 }
274
275 G4bool G4ChargeExchange::SampleMass(G4double& M, const G4double G, const G4double elim)
276 {
277 // +- 4 width but above 2 pion mass
278 const G4double e1 = std::max(M - 4*G, emin);
279 const G4double e2 = std::min(M + 4*G, elim) - e1;
280 if (e2 <= 0.0) { return false; }
281 const G4double M2 = M*M;
282 const G4double MG2 = M2*G*G;
283
284 // sampling Breit-Wigner function
285 for (G4int i=0; i<maxN; ++i) {
286 G4double e = e1 + e2*G4UniformRand();
287 G4double x = e*e - M2;
288 G4double y = MG2/(x*x + MG2);
289 if (y >= G4UniformRand()) {
290 M = e;
291 return true;
292 }
293 }
294 return false;
295 }
296