| 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 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35
36 #include "globals.hh"
37
38 /** \file G4INCLCrossSectionsAntiparticles.cc
39 * \brief Multipion, mesonic Resonances, strange cross sections and antinucleon as projectile
40 *
41 * \date 31st March 2023
42 * \author Demid Zharenov
43 */
44
45 #include "G4INCLCrossSectionsAntiparticles.hh"
46 #include "G4INCLKinematicsUtils.hh"
47 #include "G4INCLParticleTable.hh"
48 // #include <cassert>
49
50 namespace G4INCL {
51
52 template<G4int N>
53 struct BystrickyEvaluator {
54 static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients<N> const &coeffs) {
55 const G4double pMeV = pLab*1E3;
56 const G4double ekin=std::sqrt(ParticleTable::effectiveNucleonMass2+pMeV*pMeV)-ParticleTable::effectiveNucleonMass;
57 const G4double xrat=ekin*oneOverThreshold;
58 const G4double x=std::log(xrat);
59 return HornerEvaluator<N>::eval(x, coeffs) * x * std::exp(-0.5*x);
60 }
61 };
62
63 const G4int CrossSectionsAntiparticles::nMaxPiNN = 4;
64 const G4int CrossSectionsAntiparticles::nMaxPiPiN = 4;
65
66 CrossSectionsAntiparticles::CrossSectionsAntiparticles() :
67 s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
68 s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
69 s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
70 s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
71 s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
72 s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
73 s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
74 s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
75 s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
76 s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
77 {
78 }
79
80 /// \brief redefining previous cross sections
81
82 G4double CrossSectionsAntiparticles::total(Particle const * const p1, Particle const * const p2) {
83 G4double inelastic;
84 if ((p1->isNucleon() && p2->isAntiNucleon()) || (p1->isAntiNucleon() && p2->isNucleon()))
85 inelastic = NNbarCEX(p1, p2) + NNbarToNNbarpi(p1, p2) + NNbarToNNbar2pi(p1, p2) + NNbarToNNbar3pi(p1, p2) + NNbarToAnnihilation(p1, p2) + NNbarToLLbar(p1, p2);
86 else if(p1->isNucleon() && p2->isNucleon()) {
87 return CrossSectionsMultiPions::NNTot(p1, p2);
88 } else if((p1->isNucleon() && p2->isDelta()) ||
89 (p1->isDelta() && p2->isNucleon())) {
90 inelastic = CrossSectionsMultiPions::NDeltaToNN(p1, p2) + NDeltaToNLK(p1, p2) + NDeltaToNSK(p1, p2) + NDeltaToDeltaLK(p1, p2) + NDeltaToDeltaSK(p1, p2) + NDeltaToNNKKb(p1, p2);
91 } else if((p1->isNucleon() && p2->isPion()) ||
92 (p1->isPion() && p2->isNucleon())) {
93 return CrossSectionsMultiPions::piNTot(p1,p2);
94 } else if((p1->isNucleon() && p2->isEta()) ||
95 (p1->isEta() && p2->isNucleon())) {
96 inelastic = CrossSectionsMultiPionsAndResonances::etaNToPiN(p1,p2) + CrossSectionsMultiPionsAndResonances::etaNToPiPiN(p1,p2);
97 } else if((p1->isNucleon() && p2->isOmega()) ||
98 (p1->isOmega() && p2->isNucleon())) {
99 inelastic = CrossSectionsMultiPionsAndResonances::omegaNInelastic(p1,p2);
100 } else if((p1->isNucleon() && p2->isEtaPrime()) ||
101 (p1->isEtaPrime() && p2->isNucleon())) {
102 inelastic = CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(p1,p2);
103 } else if((p1->isNucleon() && p2->isLambda()) ||
104 (p1->isLambda() && p2->isNucleon())) {
105 inelastic = CrossSectionsStrangeness::NLToNS(p1,p2);
106 } else if((p1->isNucleon() && p2->isSigma()) ||
107 (p1->isSigma() && p2->isNucleon())) {
108 inelastic = CrossSectionsStrangeness::NSToNL(p1,p2) + CrossSectionsStrangeness::NSToNS(p1,p2);
109 } else if((p1->isNucleon() && p2->isKaon()) ||
110 (p1->isKaon() && p2->isNucleon())) {
111 inelastic = CrossSectionsStrangeness::NKToNK(p1,p2) + CrossSectionsStrangeness::NKToNKpi(p1,p2) + CrossSectionsStrangeness::NKToNK2pi(p1,p2);
112 } else if((p1->isNucleon() && p2->isAntiKaon()) ||
113 (p1->isAntiKaon() && p2->isNucleon())) {
114 inelastic = CrossSectionsStrangeness::NKbToLpi(p1,p2)
115 + CrossSectionsStrangeness::NKbToSpi(p1,p2) + CrossSectionsStrangeness::NKbToL2pi(p1,p2)
116 + CrossSectionsStrangeness::NKbToS2pi(p1,p2) + CrossSectionsStrangeness::NKbToNKb(p1,p2)
117 + CrossSectionsStrangeness::NKbToNKbpi(p1,p2) + CrossSectionsStrangeness::NKbToNKb2pi(p1,p2);
118 } else {
119 inelastic = 0.;
120 }
121 return inelastic + elastic(p1, p2);
122 }
123
124 // without NNbar!
125 G4double CrossSectionsAntiparticles::elastic(Particle const * const p1, Particle const * const p2) {
126 if ((p1->isNucleon() && p2->isAntiNucleon()) || (p1->isAntiNucleon() && p2->isNucleon()))
127 return NNbarElastic(p1, p2);
128 if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){ // N-N, N-Delta, Delta-Delta
129 return CrossSectionsMultiPions::elastic(p1, p2);
130 }
131 else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
132 return CrossSectionsMultiPions::elastic(p1, p2);
133 }
134 else if ((p1->isNucleon() && p2->isEta()) || (p2->isNucleon() && p1->isEta())){
135 return CrossSectionsMultiPionsAndResonances::etaNElastic(p1, p2);
136 }
137 else if ((p1->isNucleon() && p2->isHyperon()) || (p2->isNucleon() && p1->isHyperon())){
138 return CrossSectionsStrangeness::NYelastic(p1, p2);
139 }
140 else if ((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon())){
141 return CrossSectionsStrangeness::NKelastic(p1, p2);
142 }
143 else if ((p1->isNucleon() && p2->isAntiKaon()) || (p2->isNucleon() && p1->isAntiKaon())){
144 return CrossSectionsStrangeness::NKbelastic(p1, p2);
145 }
146 else {
147 return 0.0;
148 }
149 }
150
151 G4double CrossSectionsAntiparticles::NNbarCEX(Particle const * const p1, Particle const * const p2) {
152 //brief ppbar
153 // p pbar -> n nbar (BFMM 204)
154 //
155 //brief nnbar
156 // n nbar -> p pbar (same as BFMM 204, but no threshold)
157 //
158
159 // assert((p1->isAntiNucleon() && p2->isNucleon()) || (p1->isNucleon() && p2->isAntiNucleon()));
160
161 G4double sigma=0.;
162 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
163 // iso == 2 || iso == -2 (n pbar or p nbar)
164
165 const std::vector<G4double> BFMM204 = {7.549, -0.041, -2.959, -6.835, 1.629, 0.114};
166 //{6.875, 0.590, -0.003, -6.629, 1.532, 0.114}
167 //const G4double Eth_PPbar_NNbar = 0.114;
168 const std::vector<G4double> BFMM204nn = {7.549, -0.041, -2.959, -6.835, 1.629};
169 //const G4double Eth_NNbar_PPbar = 0.0;
170
171 const Particle *antinucleon;
172 const Particle *nucleon;
173
174 if (p1->isAntiNucleon()) {
175 antinucleon = p1;
176 nucleon = p2;
177 }
178 else {
179 antinucleon = p2;
180 nucleon = p1;
181 }
182
183 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antinucleon, nucleon); // GeV
184
185 if(iso == 2 || iso == -2){ //npbar or pnbar
186 sigma = 0.0;
187 return sigma;
188 }
189 else{ // ppbar or nnbar
190 if(p1->getType()==antiProton || p1->getType()==Proton)
191 sigma = KinematicsUtils::compute_xs(std::move(BFMM204), pLab); // ppbar case
192 else
193 sigma = KinematicsUtils::compute_xs(std::move(BFMM204nn), pLab); // nnbar case
194 return sigma;
195 }
196 }
197
198 G4double CrossSectionsAntiparticles::NNbarElastic(Particle const * const p1, Particle const * const p2) {
199 //brief ppbar
200 // p pbar -> p pbar (BFMM 2)
201 //
202 //brief npbar
203 // n pbar -> n pbar (BFMM 472)
204 //
205 //brief nnbar
206 // n nbar -> n nbar (same as BFMM 2)
207 //
208 //brief pnbar
209 // p nbar -> p nbar (same as BFMM 472)
210 //
211
212 // assert((p1->isAntiNucleon() && p2->isNucleon()) || (p1->isNucleon() && p2->isAntiNucleon()));
213
214 G4double sigma=0.;
215 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
216 // iso == 2 || iso == -2 (n pbar or p nbar)
217
218 const std::vector<G4double> BFMM2 = {110.496, -65.605, -0.198, -34.813, 4.317};
219 //elastic ppbar;
220 const std::vector<G4double> BFMM472 = {14.625, 23.413, -0.288, -9.002, 1.084};
221 //elastic pnbar;
222
223 const Particle *antinucleon;
224 const Particle *nucleon;
225
226 if (p1->isAntiNucleon()) {
227 antinucleon = p1;
228 nucleon = p2;
229 }
230 else {
231 antinucleon = p2;
232 nucleon = p1;
233 }
234
235 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antinucleon, nucleon); // GeV
236
237 if(iso == 2 || iso == -2){ // npbar or pnbar
238 sigma = KinematicsUtils::compute_xs(std::move(BFMM472), pLab);
239 return sigma;
240 }
241 else { // ppbar or nnbar
242 sigma = KinematicsUtils::compute_xs(std::move(BFMM2), pLab);
243 return sigma;
244 }
245 }
246
247 G4double CrossSectionsAntiparticles::NNbarToLLbar(Particle const * const p1, Particle const * const p2) {
248 // this channel includes all states with lambdas, sigmas and xis and their antiparticles
249
250 //brief ppbar
251 // p pbar -> l lbar (BFMM 121)
252 // ppbar -> l lbar pi0 (BFMM 113)
253 // ppbar -> splus pim lbar || sminusbar pim l (BFMM 136)
254 // ppbar -> sminus pip lbar || splusbar l pip (BFMM 146)
255 // ppbar -> sp spbar (BFMM 139)
256 // ppbar -> sm smbar (BFMM 149)
257 // ppbar -> szero szerobar (BFMM 144)
258 // ppbar -> ximinus ximinusbar (BFMM 101)
259 // ppbar -> szero lbar || szerobar l (BFMM 143)
260 //
261 //
262 //brief npbar
263 // n pbar -> l lbar pi- (BFMM 487)
264 // n pbar -> l sbarplus || lbar sminus (BFMM 488)
265 //
266 //
267 //brief nnbar
268 // all same as for ppbar
269 //
270 //
271 //brief pnbar
272 // p nbar -> l lbar pi+ (same as BFMM 487)
273 // p nbar -> l sbarminus || lbar splus (same as BFMM 488)
274 //
275
276 const std::vector<G4double> BFMM121 = {2.379, -2.738, -1.260, -1.915, 0.430, 1.437};
277 //const G4double Eth_PPbar_LLbar = 1.437;
278 const std::vector<G4double> BFMM113 = {-0.105, 0.000, -5.099, 0.188, -0.050, 1.820};
279 //const G4double Eth_PPbar_LLbar_pi0 = 1.820;
280 const std::vector<G4double> BFMM139 = {0.142, -0.291, -1.702, -0.058, 0.001, 1.851};
281 //const G4double Eth_PPbar_SpSpbar = 1.851;
282 const std::vector<G4double> BFMM149 = {1.855, -2.238, -1.002, -1.279, 0.252, 1.896};
283 //const G4double Eth_PPbar_SmSmbar = 1.896;
284 const std::vector<G4double> BFMM136 = {1.749, -2.506, -1.222, -1.262, 0.274, 2.042};
285 //const G4double Eth_PPbar_SpLbar_pim = 2.042;
286 const std::vector<G4double> BFMM146 = {1.037, -1.437, -1.155, -0.709, 0.138, 2.065};
287 //const G4double Eth_PPbar_SmLbar_pip = 2.065;
288 const std::vector<G4double> BFMM143 = {0.652, -1.006, -1.805, -0.537, 0.121, 1.653};
289 //const G4double Eth_PPbar_Szero_Lbar = 1.653;
290
291 // assert((p1->isAntiNucleon() && p2->isNucleon()) || (p1->isNucleon() && p2->isAntiNucleon()));
292
293 G4double sigma=0.;
294 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
295 // iso == 2 || iso == -2 (n pbar or p nbar)
296
297 const Particle *antinucleon;
298 const Particle *nucleon;
299
300 if (p1->isAntiNucleon()) {
301 antinucleon = p1;
302 nucleon = p2;
303 }
304 else {
305 antinucleon = p2;
306 nucleon = p1;
307 }
308
309 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antinucleon, nucleon); // GeV
310
311 //fixed due to limited data
312 G4double BFMM144;
313 if(pLab > 1.868) BFMM144 = 0.008; //sigmazero sigmazerobar
314 else BFMM144 = 0.0;
315 G4double BFMM101;
316 if(pLab > 1.868) BFMM101 = 0.002; //xizero xizerobar
317 else BFMM101 = 0.0;
318
319 // npbar cross sections (fixed due to limited data)
320 G4double BFMM487;
321 if(pLab > 2.1) BFMM487 = 0.048; //llbar piminus
322 else BFMM487 = 0.0;
323 G4double BFMM488;
324 if(pLab > 2.0) BFMM488 = 0.139; //lsigmaminus +cc
325 else BFMM488 = 0.0;
326
327 if(iso == 2 || iso == -2){ //npbar or pnbar
328 sigma = BFMM487 + BFMM488;
329 return sigma;
330 }
331 else{ // ppbar or nnbar
332 sigma = KinematicsUtils::compute_xs(std::move(BFMM113), pLab)
333 +KinematicsUtils::compute_xs(std::move(BFMM139), pLab)+KinematicsUtils::compute_xs(std::move(BFMM136), pLab)
334 +KinematicsUtils::compute_xs(std::move(BFMM146), pLab)+KinematicsUtils::compute_xs(std::move(BFMM143), pLab)
335 +KinematicsUtils::compute_xs(std::move(BFMM121), pLab)+KinematicsUtils::compute_xs(std::move(BFMM149), pLab)
336 +BFMM144 +BFMM101; // nnbar case totally same as ppbar
337 return sigma;
338 }
339 }
340
341 G4double CrossSectionsAntiparticles::NNbarToNNbarpi(Particle const * const p1, Particle const * const p2) {
342 //brief ppbar
343 // p pbar -> p pbar pi0 (BFMM 185)
344 // p pbar -> p nbar pi- (BFMM 188)
345 // p pbar -> n pbar pi+ (BFMM 199)
346 // p pbar -> n nbar pi0 (no data)
347 //
348 //brief npbar
349 // n pbar -> p pbar pi- (BFMM 491)
350 // n pbar -> p nbar pion (impossible)
351 // n pbar -> n pbar pi0 (BFMM 495)
352 // n pbar -> n nbar pi- (same as BFMM 188)
353 //
354 //brief nnbar
355 // n nbar -> n nbar pi0 (same as BFMM 185)
356 // n nbar -> p nbar pi- (same as BFMM 188)
357 // n nbar -> n pbar pi+ (same as BFMM 199)
358 // n nbar -> p pbar pi0 (no data)
359 //
360 //brief pnbar
361 // p nbar -> p pbar pi+ (same as BFMM 491)
362 // p nbar -> n pbar pion (impossible)
363 // p nbar -> p nbar pi0 (BFMM 495)
364 // p nbar -> n nbar pi- (same as BFMM 188)
365 //
366 //
367 // BFMM 188,199 are very close in value, 491 is larger
368
369 // assert((p1->isAntiNucleon() && p2->isNucleon()) || (p1->isNucleon() && p2->isAntiNucleon()));
370
371 G4double sigma=0.;
372 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
373 // iso == 2 || iso == -2 (n pbar or p nbar)
374
375 const std::vector<G4double> BFMM185 = {-0.734, 0.841, 0.905, 3.415, -2.316, 0.775};
376 //{22.781, -22.602, -0.752, -11.036, 1.548, 0.775}
377 //const G4double Eth_PPbar_PPbar_pi0 = 0.775;
378 const std::vector<G4double> BFMM188 = { -0.442, 0.501, 0.002, 3.434, -1.201, 0.798};
379 //const G4double Eth_PPbar_PNbar_pim = 0.798;
380 const std::vector<G4double> BFMM199 = {-2.025, 2.055, -2.355, 6.064, -2.004, 0.798};
381 //const G4double Eth_PPbar_NPbar_pip = 0.798;
382 const std::vector<G4double> BFMM491 = {24.125, -20.669, -1.534, -19.573, 4.493, 0.787};
383 //const G4double Eth_NPbar_PPbar_pim = 0.787;
384 const std::vector<G4double> BFMM495 = {-0.650, -0.140, -0.058, 5.166, -1.705, 0.777};
385 //const G4double Eth_NPbar_NPbar_pi0 = 0.777;
386
387 const Particle *antinucleon;
388 const Particle *nucleon;
389
390 if (p1->isAntiNucleon()) {
391 antinucleon = p1;
392 nucleon = p2;
393 }
394 else {
395 antinucleon = p2;
396 nucleon = p1;
397 }
398
399 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antinucleon, nucleon); // GeV
400
401 if(iso == 2 || iso == -2){ //npbar or pnbar
402 sigma = KinematicsUtils::compute_xs(std::move(BFMM491), pLab) + KinematicsUtils::compute_xs(std::move(BFMM185), pLab) + KinematicsUtils::compute_xs(std::move(BFMM188), pLab);
403 return sigma;
404 }
405 else{ // ppbar or nnbar
406 sigma = KinematicsUtils::compute_xs(std::move(BFMM199), pLab) + KinematicsUtils::compute_xs(std::move(BFMM185), pLab) + KinematicsUtils::compute_xs(std::move(BFMM188), pLab);
407 return sigma;
408 }
409 }
410
411 G4double CrossSectionsAntiparticles::NNbarToNNbar2pi(Particle const * const p1, Particle const * const p2) {
412 //brief ppbar
413 // p pbar -> p pbar pi+ pi- (BFMM 167)
414 // p pbar -> p nbar pi- pi0 (same as BFMM 490)
415 // p pbar -> n pbar pi+ pi0 (same as BFMM 490)
416 // p pbar -> n nbar pi+ pi- (BFMM 198)
417 //
418 //brief npbar
419 // n pbar -> p pbar pi- pi0 (BFMM 490)
420 // n pbar -> p nbar pi- pi- (BFMM 492)
421 // n pbar -> n pbar pi+ pi- (BFMM 494)
422 // n pbar -> n nbar pi- pi0 (same as BFMM 490)
423 //
424 //brief nnbar
425 // n nbar -> n nbar pi+ pi- (same as BFMM 167)
426 // n nbar -> p nbar pi- pi0 (same as BFMM 490)
427 // n nbar -> n pbar pi+ pi0 (same as BFMM 490)
428 // n nbar -> p pbar pi+ pi- (same as BFMM 198)
429 //
430 //brief pnbar
431 // p nbar -> p pbar pi+ pi0 (same as BFMM 490)
432 // p nbar -> n pbar pi+ pi+ (same as BFMM 492)
433 // p nbar -> p nbar pi+ pi- (same as BFMM 494)
434 // p nbar -> n nbar pi+ pi0 (same as BFMM 490)
435 //
436 //
437 // BFMM 188,199 are very close in value, 491 is larger
438
439 // assert((p1->isAntiNucleon() && p2->isNucleon()) || (p1->isNucleon() && p2->isAntiNucleon()));
440
441 G4double sigma=0.;
442 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
443 // iso == 2 || iso == -2 (n pbar or p nbar)
444
445 const std::vector<G4double> BFMM167 = {-6.885, 0.476, 1.206, 13.857, -5.728, 1.220};
446 //const G4double Eth_PPbar_PPbar_pip_pim = 1.220;
447 const std::vector<G4double> BFMM198 = {1.857, -21.213, -3.448, 0.827, -0.390, 1.231};
448 //const G4double Eth_PPbar_NNbar_pip_pim = 1.231;
449 const std::vector<G4double> BFMM490 = {-3.594, 0.811, 0.306, 5.108, -1.625, 1.201};
450 //const G4double Eth_PNbar_PPbar_pim_pi0 = 1.201;
451 const std::vector<G4double> BFMM492 = {-5.443, 7.254, -2.936, 8.441, -2.588, 1.221};
452 //const G4double Eth_PNbar_NPbar_pim_pim = 1.221;
453 const std::vector<G4double> BFMM494 = {21.688, -38.709, -2.062, -17.783, 3.895, 1.221};
454 //const G4double Eth_NPbar_NPbar_pip_pim = 1.221;
455
456 const Particle *antinucleon;
457 const Particle *nucleon;
458
459 if (p1->isAntiNucleon()) {
460 antinucleon = p1;
461 nucleon = p2;
462 }
463 else {
464 antinucleon = p2;
465 nucleon = p1;
466 }
467
468 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antinucleon, nucleon); // GeV
469
470 if(iso == 2 || iso == -2){ // pnbar or npbar
471 sigma = KinematicsUtils::compute_xs(BFMM490, pLab) + KinematicsUtils::compute_xs(BFMM490, pLab) + KinematicsUtils::compute_xs(std::move(BFMM167), pLab) + KinematicsUtils::compute_xs(std::move(BFMM198), pLab);
472 return sigma;
473 }
474 else{ // ppbar or nnbar
475 sigma = KinematicsUtils::compute_xs(BFMM490, pLab) + KinematicsUtils::compute_xs(BFMM490, pLab) + KinematicsUtils::compute_xs(std::move(BFMM492), pLab) + KinematicsUtils::compute_xs(std::move(BFMM494), pLab);
476 return sigma;
477 }
478 }
479
480 G4double CrossSectionsAntiparticles::NNbarToNNbar3pi(Particle const * const p1, Particle const * const p2) {
481 //brief ppbar
482 // p pbar -> p pbar pi+ pi- pi0 (BFMM 161)
483 // p pbar -> p nbar 2pi- pi+ (BFMM 169)
484 // p pbar -> n pbar 2pi+ pi- (BFMM 201)
485 // p pbar -> n nbar pi+ pi- pi0 (BFMM 197)
486 //
487 //brief npbar
488 // n pbar -> p pbar 2pi- pi+ (same as BFMM 169)
489 // n pbar -> p nbar 2pi- pi0 (same as BFMM 197)
490 // n pbar -> n pbar pi+ pi- pi0 (same as BFMM 161)
491 // n pbar -> n nbar 2pi- pi+ (same as BFMM 169)
492 //
493 //brief nnbar
494 // n nbar -> n nbar pi+ pi- pi0 (same as BFMM 161)
495 // n nbar -> p nbar 2pi- pi+ (same as BFMM 169)
496 // n nbar -> n pbar 2pi+ pi- (same as BFMM 201)
497 // n nbar -> p pbar pi+ pi- pi0 (same as BFMM 197)
498 //
499 //brief pnbar
500 // p nbar -> p pbar 2pi+ pi- (same as BFMM 169)
501 // p nbar -> n pbar 2pi+ pi0 (same as BFMM 197)
502 // p nbar -> p nbar pi+ pi- pi0 (same as BFMM 161)
503 // p nbar -> n nbar 2pi+ pi- (same as BFMM 169)
504 //
505
506 // assert((p1->isAntiNucleon() && p2->isNucleon()) || (p1->isNucleon() && p2->isAntiNucleon()));
507
508 G4double sigma=0.;
509 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
510 // iso == 2 || iso == -2 (n pbar or p nbar)
511
512 const std::vector<G4double> BFMM161 = {-6.434, 1.351, -5.185, 7.754, -1.692, 1.604};
513 //const G4double Eth_PPbar_PPbar_pip_pim_pi0 = 1.604;
514 const std::vector<G4double> BFMM169 = {3.696, -5.356, -0.053, 1.941, -0.432, 1.624};
515 //const G4double Eth_PPbar_PNbar_2pim_pip = 1.624;
516 const std::vector<G4double> BFMM201 = {-1.070, -0.636, -0.009, 2.335, -0.499, 1.624};
517 //const G4double Eth_PPbar_NPbar_2pip_pim = 1.624;
518 const std::vector<G4double> BFMM197 = {1.857, -21.213, -3.448, 0.827, -0.390, 1.616};
519 //const G4double Eth_PPbar_NNbar_pip_pim_pi0 = 1.616;
520
521 const Particle *antinucleon;
522 const Particle *nucleon;
523
524 if (p1->isAntiNucleon()) {
525 antinucleon = p1;
526 nucleon = p2;
527 }
528 else {
529 antinucleon = p2;
530 nucleon = p1;
531 }
532
533 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antinucleon, nucleon); // GeV
534
535 if(iso == 2 || iso == -2){ // pnbar or npbar
536 sigma = KinematicsUtils::compute_xs(BFMM169, pLab) + KinematicsUtils::compute_xs(BFMM169, pLab) + KinematicsUtils::compute_xs(std::move(BFMM197), pLab) + KinematicsUtils::compute_xs(std::move(BFMM161), pLab);
537 return sigma;
538 }
539 else{ // ppbar or nnbar
540 sigma = KinematicsUtils::compute_xs(std::move(BFMM161), pLab) + KinematicsUtils::compute_xs(std::move(BFMM169), pLab) + KinematicsUtils::compute_xs(std::move(BFMM197), pLab) + KinematicsUtils::compute_xs(std::move(BFMM201), pLab);
541 return sigma;
542 }
543 }
544
545 G4double CrossSectionsAntiparticles::NNbarToAnnihilation(Particle const * const p1, Particle const * const p2) {
546 //brief ppbar
547 /*
548 This part only contains total annihilation xs, the choice of a particular final state
549 will be done in the channel file.
550 As long as we only have good data for ppbar, we assume that for npbar, pnbar and nnbar the xs
551 will be the same, but in order to compensate for the Coulombic effect the ppbar annihilation xs
552 is multiplied by the pnbar total xs and divided by the ppbar total xs.
553 */
554
555 // assert((p1->isAntiNucleon() && p2->isNucleon()) || (p1->isNucleon() && p2->isAntiNucleon()));
556
557 G4double sigma=0.;
558 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
559 // iso == 2 || iso == -2 (n pbar or p nbar)
560
561 const std::vector<G4double> BFMM6 = {66.098, 0.153, -4.576, -38.319, 6.625}; //ppbar annihilation xs
562 const std::vector<G4double> BFMM1 = {119.066, 6.251, -0.006, -60.046, 11.958}; //ppbar total xs
563 const std::vector<G4double> BFMM471 = {108.104, 15.708, 0.832, -54.632, -6.958}; //npbar total xs
564
565 const Particle *antinucleon;
566 const Particle *nucleon;
567
568 if (p1->isAntiNucleon()) {
569 antinucleon = p1;
570 nucleon = p2;
571 }
572 else {
573 antinucleon = p2;
574 nucleon = p1;
575 }
576
577 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antinucleon, nucleon); // GeV
578
579 if(iso == 2 || iso == -2){ // pnbar or npbar
580 sigma = KinematicsUtils::compute_xs(std::move(BFMM6), pLab)*KinematicsUtils::compute_xs(std::move(BFMM471), pLab)/KinematicsUtils::compute_xs(std::move(BFMM1), pLab);
581 return sigma;
582 }
583 else if(p1->getType()==antiProton || p2->getType()==Proton){ // ppbar case
584 sigma = KinematicsUtils::compute_xs(std::move(BFMM6), pLab);
585 return sigma;
586 }
587 else{ // nnbar case
588 sigma = KinematicsUtils::compute_xs(std::move(BFMM6), pLab)*KinematicsUtils::compute_xs(std::move(BFMM471), pLab)/KinematicsUtils::compute_xs(std::move(BFMM1), pLab);
589 return sigma;
590 }
591 }
592
593 } // namespace G4INCL
594
595