| 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 G4INCLCrossSectionsStrangeness.cc
39 * \brief Multipion, mesonic Resonances and strange cross sections
40 *
41 * \date 1st March 2016
42 * \author Jason Hirtz
43 */
44
45 #include "G4INCLCrossSectionsStrangeness.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 CrossSectionsStrangeness::nMaxPiNN = 4;
64 const G4int CrossSectionsStrangeness::nMaxPiPiN = 4;
65
66 CrossSectionsStrangeness::CrossSectionsStrangeness() :
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 CrossSectionsStrangeness::total(Particle const * const p1, Particle const * const p2) {
83 G4double inelastic;
84 if(p1->isNucleon() && p2->isNucleon()) {
85 return CrossSectionsMultiPions::NNTot(p1, p2);
86 } else if((p1->isNucleon() && p2->isDelta()) ||
87 (p1->isDelta() && p2->isNucleon())) {
88 inelastic = CrossSectionsMultiPions::NDeltaToNN(p1, p2) + NDeltaToNLK(p1, p2) + NDeltaToNSK(p1, p2) + NDeltaToDeltaLK(p1, p2) + NDeltaToDeltaSK(p1, p2) + NDeltaToNNKKb(p1, p2);
89 } else if((p1->isNucleon() && p2->isPion()) ||
90 (p1->isPion() && p2->isNucleon())) {
91 return CrossSectionsMultiPions::piNTot(p1,p2);
92 } else if((p1->isNucleon() && p2->isEta()) ||
93 (p1->isEta() && p2->isNucleon())) {
94 inelastic = CrossSectionsMultiPionsAndResonances::etaNToPiN(p1,p2) + CrossSectionsMultiPionsAndResonances::etaNToPiPiN(p1,p2);
95 } else if((p1->isNucleon() && p2->isOmega()) ||
96 (p1->isOmega() && p2->isNucleon())) {
97 inelastic = CrossSectionsMultiPionsAndResonances::omegaNInelastic(p1,p2);
98 } else if((p1->isNucleon() && p2->isEtaPrime()) ||
99 (p1->isEtaPrime() && p2->isNucleon())) {
100 inelastic = CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(p1,p2);
101 } else if((p1->isNucleon() && p2->isLambda()) ||
102 (p1->isLambda() && p2->isNucleon())) {
103 inelastic = NLToNS(p1,p2);
104 } else if((p1->isNucleon() && p2->isSigma()) ||
105 (p1->isSigma() && p2->isNucleon())) {
106 inelastic = NSToNL(p1,p2) + NSToNS(p1,p2);
107 } else if((p1->isNucleon() && p2->isKaon()) ||
108 (p1->isKaon() && p2->isNucleon())) {
109 inelastic = NKToNK(p1,p2) + NKToNKpi(p1,p2) + NKToNK2pi(p1,p2);
110 } else if((p1->isNucleon() && p2->isAntiKaon()) ||
111 (p1->isAntiKaon() && p2->isNucleon())) {
112 inelastic = NKbToLpi(p1,p2) + NKbToSpi(p1,p2) + NKbToL2pi(p1,p2) + NKbToS2pi(p1,p2) + NKbToNKb(p1,p2) + NKbToNKbpi(p1,p2) + NKbToNKb2pi(p1,p2);
113 } else {
114 inelastic = 0.;
115 }
116 return inelastic + elastic(p1, p2);
117 }
118
119 G4double CrossSectionsStrangeness::elastic(Particle const * const p1, Particle const * const p2) {
120 if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){ // N-N, N-Delta, Delta-Delta
121 return CrossSectionsMultiPions::elastic(p1, p2);
122 }
123 else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
124 return CrossSectionsMultiPions::elastic(p1, p2);
125 }
126 else if ((p1->isNucleon() && p2->isEta()) || (p2->isNucleon() && p1->isEta())){
127 return CrossSectionsMultiPionsAndResonances::etaNElastic(p1, p2);
128 }
129 else if ((p1->isNucleon() && p2->isHyperon()) || (p2->isNucleon() && p1->isHyperon())){
130 return NYelastic(p1, p2);
131 }
132 else if ((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon())){
133 return NKelastic(p1, p2);
134 }
135 else if ((p1->isNucleon() && p2->isAntiKaon()) || (p2->isNucleon() && p1->isAntiKaon())){
136 return NKbelastic(p1, p2);
137 }
138 else {
139 return 0.0;
140 }
141 }
142
143 G4double CrossSectionsStrangeness::piNToxPiN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
144 //
145 // pion-Nucleon producing xpi pions cross sections (corrected due to new particles)
146 //
147 // assert(xpi>1 && xpi<=nMaxPiPiN);
148 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
149
150 const G4double oldXS2Pi=CrossSectionsMultiPions::piNToxPiN(2,particle1, particle2);
151 const G4double oldXS3Pi=CrossSectionsMultiPions::piNToxPiN(3,particle1, particle2);
152 const G4double oldXS4Pi=CrossSectionsMultiPions::piNToxPiN(4,particle1, particle2);
153 const G4double xsEta=CrossSectionsMultiPionsAndResonances::piNToEtaN(particle1, particle2);
154 const G4double xsOmega=CrossSectionsMultiPionsAndResonances::piNToOmegaN(particle1, particle2);
155 const G4double xs1=NpiToLK(particle2, particle1);
156 const G4double xs2=NpiToSK(particle1, particle2);
157 const G4double xs3=NpiToLKpi(particle1, particle2);
158 const G4double xs4=NpiToSKpi(particle1, particle2);
159 const G4double xs5=NpiToLK2pi(particle1, particle2);
160 const G4double xs6=NpiToSK2pi(particle1, particle2);
161 const G4double xs7=NpiToNKKb(particle1, particle2);
162 const G4double xs8=NpiToMissingStrangeness(particle1, particle2);
163 const G4double xs0 = xs1 + xs2 + xs3 + xs4 + xs5 + xs6 + xs7 +xs8;
164 G4double newXS2Pi=0.;
165 G4double newXS3Pi=0.;
166 G4double newXS4Pi=0.;
167
168 if (xpi == 2) {
169 if (oldXS4Pi != 0.)
170 newXS2Pi=oldXS2Pi;
171 else if (oldXS3Pi != 0.) {
172 newXS3Pi=oldXS3Pi-xsEta-xsOmega-xs0;
173 if (newXS3Pi < 1.e-09)
174 newXS2Pi=oldXS2Pi-(xsEta+xsOmega+xs0-oldXS3Pi);
175 else
176 newXS2Pi=oldXS2Pi;
177 }
178 else {
179 newXS2Pi=oldXS2Pi-xsEta-xsOmega-xs0;
180 if (newXS2Pi < 1.e-09 && newXS2Pi!=0.){
181 newXS2Pi=0.;
182 }
183 }
184 return newXS2Pi;
185 }
186 else if (xpi == 3) {
187 if (oldXS4Pi != 0.) {
188 newXS4Pi=oldXS4Pi-xsEta-xsOmega-xs0;
189 if (newXS4Pi < 1.e-09)
190 newXS3Pi=oldXS3Pi-(xsEta+xsOmega+xs0-oldXS4Pi);
191 else
192 newXS3Pi=oldXS3Pi;
193 }
194 else {
195 newXS3Pi=oldXS3Pi-xsEta-xsOmega-xs0;
196 if (newXS3Pi < 1.e-09){
197 newXS3Pi=0.;
198 }
199 }
200 return newXS3Pi;
201 }
202 else if (xpi == 4) {
203 newXS4Pi=oldXS4Pi-xsEta-xsOmega-xs0;
204 if (newXS4Pi < 1.e-09){
205 newXS4Pi=0.;
206 }
207 return newXS4Pi;
208 }
209 else // should never reach this point
210 return 0.;
211 }
212
213 G4double CrossSectionsStrangeness::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
214 //
215 // Nucleon-Nucleon producing xpi pions cross sections corrected
216 //
217 // assert(xpi>0 && xpi<=nMaxPiNN);
218 // assert(particle1->isNucleon() && particle2->isNucleon());
219
220 G4double oldXS1Pi=CrossSectionsMultiPions::NNToxPiNN(1,particle1, particle2);
221 G4double oldXS2Pi=CrossSectionsMultiPions::NNToxPiNN(2,particle1, particle2);
222 G4double oldXS3Pi=CrossSectionsMultiPions::NNToxPiNN(3,particle1, particle2);
223 G4double oldXS4Pi=CrossSectionsMultiPions::NNToxPiNN(4,particle1, particle2);
224 G4double xsEta=CrossSectionsMultiPionsAndResonances::NNToNNEta(particle1, particle2)+CrossSectionsMultiPionsAndResonances::NNToNNOmega(particle1, particle2);;
225
226 const G4double xs1=NNToNLK(particle1, particle2);
227 const G4double xs2=NNToNSK(particle1, particle2);
228 const G4double xs3=NNToNLKpi(particle1, particle2);
229 const G4double xs4=NNToNSKpi(particle1, particle2);
230 const G4double xs5=NNToNLK2pi(particle1, particle2);
231 const G4double xs6=NNToNSK2pi(particle1, particle2);
232 const G4double xs7=NNToNNKKb(particle1, particle2);
233 const G4double xs8=NNToMissingStrangeness(particle1, particle2);
234 const G4double xs0 = xs1 + xs2 + xs3 + xs4 + xs5 + xs6 + xs7 + xs8;
235 G4double newXS1Pi=0.;
236 G4double newXS2Pi=0.;
237 G4double newXS3Pi=0.;
238 G4double newXS4Pi=0.;
239
240 if (xpi == 1) {
241 if (oldXS4Pi != 0. || oldXS3Pi != 0.)
242 newXS1Pi=oldXS1Pi;
243 else if (oldXS2Pi != 0.) {
244 newXS2Pi=oldXS2Pi-xsEta-xs0;
245 if (newXS2Pi < 0.)
246 newXS1Pi=oldXS1Pi-(xsEta+xs0-oldXS2Pi);
247 else
248 newXS1Pi=oldXS1Pi;
249 }
250 else
251 newXS1Pi=oldXS1Pi-xsEta-xs0;
252 return newXS1Pi;
253 }
254 else if (xpi == 2) {
255 if (oldXS4Pi != 0.){
256 newXS2Pi=oldXS2Pi;
257 }
258 else if (oldXS3Pi != 0.) {
259 newXS3Pi=oldXS3Pi-xsEta-xs0;
260 if (newXS3Pi < 0.)
261 newXS2Pi = oldXS2Pi-(xsEta+xs0-oldXS3Pi);
262 else
263 newXS2Pi = oldXS2Pi;
264 }
265 else {
266 newXS2Pi = oldXS2Pi-xsEta-xs0;
267 if (newXS2Pi < 0.)
268 newXS2Pi = 0.;
269 }
270 return newXS2Pi;
271 }
272 else if (xpi == 3) {
273 if (oldXS4Pi != 0.) {
274 newXS4Pi=oldXS4Pi-xsEta-xs0;
275 if (newXS4Pi < 0.)
276 newXS3Pi=oldXS3Pi-(xsEta+xs0-oldXS4Pi);
277 else
278 newXS3Pi=oldXS3Pi;
279 }
280 else {
281 newXS3Pi=oldXS3Pi-xsEta-xs0;
282 if (newXS3Pi < 0.)
283 newXS3Pi=0.;
284 }
285 return newXS3Pi;
286 }
287 else if (xpi == 4) {
288 newXS4Pi=oldXS4Pi-xsEta-xs0;
289 if (newXS4Pi < 0.)
290 newXS4Pi=0.;
291 return newXS4Pi;
292 }
293
294 else // should never reach this point
295 return 0.;
296 }
297
298 /// \brief elastic cross sections
299
300 G4double CrossSectionsStrangeness::NYelastic(Particle const * const p1, Particle const * const p2) {
301 //
302 // Hyperon-Nucleon elastic cross sections
303 //
304 // assert((p1->isNucleon() && p2->isHyperon()) || (p1->isHyperon() && p2->isNucleon()));
305
306 G4double sigma = 0.;
307
308 const Particle *hyperon;
309 const Particle *nucleon;
310
311 if (p1->isHyperon()) {
312 hyperon = p1;
313 nucleon = p2;
314 }
315 else {
316 hyperon = p2;
317 nucleon = p1;
318 }
319
320 const G4double pLab = KinematicsUtils::momentumInLab(hyperon, nucleon); // MeV
321
322 if (pLab < 145.)
323 sigma = 200.;
324 else if (pLab < 425.)
325 sigma = 869.*std::exp(-pLab/100.);
326 else if (pLab < 30000.)
327 sigma = 12.8*std::exp(-6.2e-5*pLab);
328 else
329 sigma=0.;
330
331 if (sigma < 0.) sigma = 0.; // should never happen
332 return sigma;
333 }
334
335 G4double CrossSectionsStrangeness::NKelastic(Particle const * const p1, Particle const * const p2) {
336 //
337 // Kaon-Nucleon elastic cross sections
338 //
339 // assert((p1->isNucleon() && p2->isKaon()) || (p1->isKaon() && p2->isNucleon()));
340
341 G4double sigma=0.;
342
343 const Particle *kaon;
344 const Particle *nucleon;
345
346 if (p1->isKaon()) {
347 kaon = p1;
348 nucleon = p2;
349 }
350 else {
351 kaon = p2;
352 nucleon = p1;
353 }
354
355 const G4double pLab = KinematicsUtils::momentumInLab(kaon, nucleon); // MeV
356
357 if (pLab < 935.)
358 sigma = 12.;
359 else if (pLab < 2080.)
360 sigma = 17.4-3.*std::exp(6.3e-4*pLab);
361 else if (pLab < 5500.)
362 sigma = 832.*std::pow(pLab,-0.64);
363 else if (pLab < 30000.)
364 sigma = 3.36;
365 else
366 sigma=0.;
367
368 if (sigma < 0.) sigma = 0.; // should never happen
369 return sigma;
370 }
371
372 G4double CrossSectionsStrangeness::NKbelastic(Particle const * const p1, Particle const * const p2) {
373 //
374 // antiKaon-Nucleon elastic cross sections
375 //
376 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
377
378 G4double sigma=0.;
379
380 const Particle *antikaon;
381 const Particle *nucleon;
382
383 if (p1->isAntiKaon()) {
384 antikaon = p1;
385 nucleon = p2;
386 }
387 else {
388 antikaon = p2;
389 nucleon = p1;
390 }
391
392 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
393
394 if(pLab > 1E-6) // sigma = 287.823 [mb] -> rise very slowly, not cut needed
395 sigma = 6.132*std::pow(pLab,-0.2437)+12.98*std::exp(-std::pow(pLab-0.9902,2)/0.05558)+2.928*std::exp(-std::pow(pLab-1.649,2)/0.772)+564.3*std::exp(-std::pow(pLab+0.9901,2)/0.5995);
396
397 if (sigma < 0.) sigma = 0.; // should never happen
398 return sigma;
399 }
400
401 /// \brief NN to strange cross sections
402
403 G4double CrossSectionsStrangeness::NNToNLK(Particle const * const p1, Particle const * const p2) {
404 //
405 // Nucleon-Nucleon producing N-Lambda-Kaon cross sections
406 //
407 // ratio
408 // p p (1) p n (1)
409 //
410 // p p -> p L K+ (1)
411 // p n -> p L K0 (1/2)
412 // p n -> n L K+ (1/2)
413 // assert(p1->isNucleon() && p2->isNucleon());
414
415 const Particle *particle1;
416 const Particle *particle2;
417
418 if(p2->getType() == Proton && p1->getType() == Neutron){
419 particle1 = p2;
420 particle2 = p1;
421 }
422 else{
423 particle1 = p1;
424 particle2 = p2;
425 }
426
427 G4double sigma = 0.;
428
429 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1, particle2); // GeV
430
431 if(particle2->getType() == Proton){
432 if(pLab < 2.3393) return 0.;
433 else if (pLab < 30.) sigma = 1.11875*std::pow((pLab-2.3393),1.0951)/std::pow((pLab+2.3393),2.0958); // pLab = 30 Gev -> exess of energie = 5 GeV
434 else return 0.;
435 }
436 else{
437 if(pLab < 2.3508) return 0.;
438 else if (pLab < 30.) sigma = 1.11875*std::pow((pLab-2.3508),1.0951)/std::pow((pLab+2.3508),2.0958);
439 else return 0.;
440 }
441
442 return sigma;
443 }
444
445 G4double CrossSectionsStrangeness::NNToNSK(Particle const * const p1, Particle const * const p2) {
446 //
447 // Nucleon-Nucleon producing N-Sigma-Kaon cross sections
448 //
449 // Meson symmetry
450 // pp->pS+K0 (1/4)
451 // pp->pS0K+ (1/8) // HEM
452 // pp->pS0K+ (1/4) // Data
453 // pp->nS+K+ (1)
454 //
455 // pn->nS+K0 (1/4)
456 // pn->pS-K+ (1/4)
457 // pn->nS0K+ (5/8)
458 // pn->pS0K0 (5/8)
459 //
460 // assert(p1->isNucleon() && p2->isNucleon());
461
462 const Particle *particle1;
463 const Particle *particle2;
464
465 if(p2->getType() == Proton && p1->getType() == Neutron){
466 particle1 = p2;
467 particle2 = p1;
468 }
469 else{
470 particle1 = p1;
471 particle2 = p2;
472 }
473
474 G4double sigma = 0.;
475
476 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1, particle2); // GeV
477
478 if(pLab < 2.593)
479 return 0.;
480
481 if(p2->getType() == p1->getType())
482 // sigma = 1.375*2*6.38*std::pow(pLab-2.57,2.1)/std::pow(pLab,4.162);
483 sigma = 1.5*6.38*std::pow(pLab-2.593,2.1)/std::pow(pLab,4.162);
484 else
485 sigma = 1.75*6.38*std::pow(pLab-2.593,2.1)/std::pow(pLab,4.162);
486
487
488 return sigma;
489 }
490
491 G4double CrossSectionsStrangeness::NNToNLKpi(Particle const * const p1, Particle const * const p2) {
492 //
493 // Nucleon-Nucleon producing N-Lambda-Kaon-pion cross sections
494 //
495 // ratio (pure NN -> DLK)
496 // pp (12) pn (8)
497 //
498 // pp -> p pi+ L K0 (9)(3)
499 // pp -> p pi0 L K+ (2)(1*2/3)
500 // pp -> n pi+ L K+ (1)(1*1/3)
501 //
502 // pn -> p pi- L K+ (2)(2*1/3)
503 // pn -> n pi0 L K+ (4)(2*2/3)
504 // pn -> p pi0 L K0 (4)
505 // pn -> n pi+ L K0 (2)
506
507 // assert(p1->isNucleon() && p2->isNucleon());
508
509 G4double sigma = 0.;
510 G4double ratio = 0.;
511 G4double ratio1 = 0.;
512 G4double ratio2 = 0.;
513 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 540.;
514 if( ener < p1->getMass() + p2->getMass())
515 return 0;
516 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
517
518 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
519 if (iso != 0){
520 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
521 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
522 }
523 else {
524 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
525 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
526 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
527 }
528
529 if( ratio1 == 0 || ratio2 == 0)
530 return 0.;
531
532 ratio = ratio2/ratio1;
533
534 sigma = ratio * NNToNLK(p1,p2) * 3;
535
536 /* const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1, p2); // GeV
537 if(pLab <= 2.77) return 0.;
538 sigma = 0.4 * std::pow(pLab-2.77,1.603)/std::pow(pLab,1.492);*/
539
540 return sigma;
541 }
542
543 G4double CrossSectionsStrangeness::NNToNSKpi(Particle const * const p1, Particle const * const p2) {
544 //
545 // Nucleon-Nucleon producing N-Sigma-Kaon-pion cross sections
546 //
547 // ratio (pure NN -> DSK)
548 // pp (36) pn (36)
549 //
550 // pp -> p pi+ S- K+ (9)
551 // pp -> p pi+ S0 K0 (9)
552 // pp -> p pi0 S+ K0 (4)
553 // pp -> n pi+ S+ K0 (2)
554 // pp -> p pi0 S0 K+ (4)
555 // pp -> n pi+ S0 K+ (2)
556 // pp -> p pi- S+ K+ (2)
557 // pp -> n pi0 S+ K+ (4)
558
559 // pn -> p pi0 S- K+ (4)
560 // pn -> n pi+ S- K+ (2)
561 // pn -> p pi0 S0 K0 (2)
562 // pn -> n pi+ S0 K0 (1)
563 // pn -> p pi+ S- K0 (9)
564
565 // assert(p1->isNucleon() && p2->isNucleon());
566
567 G4double sigma = 0.;
568 G4double ratio = 0.;
569 G4double ratio1 = 0.;
570 G4double ratio2 = 0.;
571 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 620.;
572 if( ener < p1->getMass() + p2->getMass())
573 return 0;
574 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
575
576 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
577 if (iso != 0){
578 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
579 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
580 }
581 else {
582 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
583 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
584 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
585 }
586
587 if( ratio1 == 0 || ratio2 == 0)
588 return 0.;
589
590 ratio = ratio2/ratio1;
591
592 sigma = ratio * NNToNSK(p1,p2) * 3;
593
594 return sigma;
595 }
596
597 G4double CrossSectionsStrangeness::NNToNLK2pi(Particle const * const p1, Particle const * const p2) {
598 //
599 // Nucleon-Nucleon producing N-Lambda-Kaon-pion cross sections
600 //
601 // assert(p1->isNucleon() && p2->isNucleon());
602
603 G4double sigma = 0.;
604 G4double ratio = 0.;
605 G4double ratio1 = 0.;
606 G4double ratio2 = 0.;
607 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 675.;
608 if( ener < p1->getMass() + p2->getMass())
609 return 0;
610 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
611
612 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
613 if (iso != 0){
614 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
615 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
616 }
617 else {
618 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
619 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
620 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
621 }
622
623 if( ratio1 == 0 || ratio2 == 0)
624 return 0.;
625
626 ratio = ratio2/ratio1;
627
628 sigma = ratio * NNToNLKpi(p1,p2);
629
630 return sigma;
631 }
632
633 G4double CrossSectionsStrangeness::NNToNSK2pi(Particle const * const p1, Particle const * const p2) {
634 //
635 // Nucleon-Nucleon producing N-Sigma-Kaon-pion cross sections
636 //
637 // assert(p1->isNucleon() && p2->isNucleon());
638
639 G4double sigma = 0.;
640 G4double ratio = 0.;
641 G4double ratio1 = 0.;
642 G4double ratio2 = 0.;
643 const G4double ener = KinematicsUtils::totalEnergyInCM(p1, p2) - 755.;
644 if( ener < p1->getMass() + p2->getMass())
645 return 0;
646 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
647
648 const G4double xsiso2 = CrossSectionsMultiPions::NNInelasticIso(ener, 2);
649 if (iso != 0){
650 ratio1 = CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
651 ratio2 = CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
652 }
653 else {
654 const G4double xsiso0 = CrossSectionsMultiPions::NNInelasticIso(ener, 0);
655 ratio1 = 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
656 ratio2 = 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
657 }
658
659 if( ratio1 == 0 || ratio2 == 0)
660 return 0.;
661
662 ratio = ratio2/ratio1;
663
664 sigma = ratio * NNToNSKpi(p1,p2);
665
666 return sigma;
667 }
668
669 G4double CrossSectionsStrangeness::NNToNNKKb(Particle const * const p1, Particle const * const p2) {
670 //
671 // Nucleon-Nucleon producing Nucleon-Nucleon-Kaon-antiKaon cross sections
672 //
673 // Channel strongly resonant; fit from Sibirtesev - Z. Phys. A 358, 101-106 (1997) (eq.21)
674 // ratio
675 // pp (6) pn (13)*2
676 // pp -> pp K+ K- (1)
677 // pp -> pp K0 K0 (1)
678 // pp -> pn K+ K0 (4)
679 // pn -> pp K0 K- (4)
680 // pn -> pn K+ K- (9)
681 //
682
683 // assert(p1->isNucleon() && p2->isNucleon());
684
685 G4double sigma = 0.;
686 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
687 const G4double ener = 0.001*KinematicsUtils::totalEnergyInCM(p1, p2); // GeV
688
689 if(ener < 2.872)
690 return 0.;
691
692 if(iso == 0)
693 sigma = 26 * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
694 else
695 sigma = 6 * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
696
697 return sigma;
698 }
699
700 G4double CrossSectionsStrangeness::NNToMissingStrangeness(Particle const * const p1, Particle const * const p2) {
701 //
702 // Nucleon-Nucleon missing strangeness production cross sections
703 //
704 // assert(p1->isNucleon() && p2->isNucleon());
705
706 G4double sigma = 0.;
707
708 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
709 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
710
711 if(pLab < 6.) return 0.;
712
713 if(iso == 0){
714 if(pLab < 30.) sigma = 10.15*std::pow((pLab - 6.),2.157)/std::pow(pLab,2.333);
715 else return 0.;
716 }
717 else{
718 if(pLab < 30.) sigma = 8.12*std::pow((pLab - 6.),2.157)/std::pow(pLab,2.333);
719 else return 0.;
720 }
721 return sigma;
722 }
723
724 /** \brief NDelta to strange cross sections
725 *
726 * No experimental data
727 * Parametrization from Phys.Rev.C 59 1 (369) (1999)
728 *
729 * Correction are applied on the isospin symetry provided in the paper
730 */
731
732 G4double CrossSectionsStrangeness::NDeltaToNLK(Particle const * const p1, Particle const * const p2) {
733 // Nucleon-Delta producing Nucleon Lambda Kaon cross section
734 //
735 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
736 //
737 // ratio
738 // D++ n -> p L K+ (3)
739 //
740 // D+ p -> p L K+ (1)
741 //
742 // D+ n -> p L K0 (1)
743 // D+ n -> n L K+ (1)
744
745 G4double a = 4.169;
746 G4double b = 2.227;
747 G4double c = 2.511;
748 G4double n_channel = 4.; // number of channel divided by 2. Here 8/2
749
750 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
751
752 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
753 if(std::abs(iso) == 4) return 0.;
754
755 G4double sigma = 0.;
756
757 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // MeV^2
758 const G4double s0 = 6.511E6; // MeV^2
759
760 if(s <= s0) return 0.;
761
762 sigma = n_channel*a*std::pow(s/s0-1,b)*std::pow(s0/s,c);
763
764 //const G4double pLab = sdt::sqrt(s*s/(4*ParticleTable::effectiveNucleonMass2)-s)*0.001;
765 //sigma = 3*1.11875*std::pow((pLab-2.3508),1.0951)/std::pow((pLab+2.3508),2.0958); // NDelta sim to NN
766
767 if(iso == 0){ // D+ n
768 sigma *= 2./6.;
769 }
770 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType())){// D+ p
771 sigma *= 1./6.;
772 }
773 else{// D++ n
774 sigma *= 3./6.;
775 }
776 return sigma;
777 }
778 G4double CrossSectionsStrangeness::NDeltaToNSK(Particle const * const p1, Particle const * const p2) {
779 // Nucleon-Delta producing Nucleon Sigma Kaon cross section
780 //
781 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369 ( X 1.25 (124/99) for isospin consideration)
782 //
783 // ratio
784 // D++ p -> p S+ K+ (6)
785 //
786 // D++ n -> p S+ K0 (3) ****
787 // D++ n -> p S0 K+ (3)
788 // D++ n -> n S+ K+ (3)
789 //
790 // D+ p -> p S+ K0 (2)
791 // D+ p -> p S0 K+ (2)
792 // D+ p -> n S+ K+ (3)
793 //
794 // D+ n -> p S0 K0 (3)
795 // D+ n -> p S- K+ (2)
796 // D+ n -> n S+ K0 (2)
797 // D+ n -> n S0 K+ (2)
798
799 G4double a = 39.54;
800 G4double b = 2.799;
801 G4double c = 6.303;
802 G4double n_channel = 11.;
803
804 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
805
806 G4double sigma = 0.;
807
808 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
809 const G4double s0 = 6.935E6; // Mev^2
810 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
811
812 if(s <= s0)
813 return 0.;
814
815 sigma = n_channel*a*std::pow(s/s0-1,b)*std::pow(s0/s,c);
816
817 //const G4double pLab = sdt::sqrt(s*s/(4*ParticleTable::effectiveNucleonMass2)-s)*0.001;
818 //sigma = 22./12./2. * 4.75*6.38*std::pow(pLab-2.593,2.1)/std::pow(pLab,4.162); // NDelta sim to NN
819
820 if(iso == 0)// D+ n
821 sigma *= 9./31.;
822 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
823 sigma *= 7./31.;
824 else if (std::abs(iso) == 2)// D++ n
825 sigma *= 9./31.;
826 else // D++ p
827 sigma *= 6./31.;
828
829 return sigma;
830 }
831 G4double CrossSectionsStrangeness::NDeltaToDeltaLK(Particle const * const p1, Particle const * const p2) {
832 // Nucleon-Delta producing Delta Lambda Kaon cross section
833 //
834 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
835 //
836 // ratio
837 // D++ p -> L K+ D++ (4)
838 //
839 // D++ n -> L K+ D+ (3)
840 // D++ n -> L K0 D++ (4)
841 //
842 // D+ p -> L K0 D++ (3)
843 // D+ p -> L K+ D+ (2)
844 //
845 // D+ n -> L K+ D0 (4)
846 // D+ n -> L K0 D+ (2)
847
848 G4double a = 2.679;
849 G4double b = 2.280;
850 G4double c = 5.086;
851 G4double n_channel = 7.;
852
853 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
854
855 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
856 const G4double s0 = 8.096E6; // Mev^2
857 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
858
859 if(s <= s0)
860 return 0.;
861
862 G4double sigma = n_channel*a*std::pow(s/s0-1,b)*std::pow(s0/s,c);
863
864 if(iso == 0)// D+ n
865 sigma *= 6./22.;
866 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
867 sigma *= 5./22.;
868 else if (std::abs(iso) == 2)// D++ n
869 sigma *= 7./22.;
870 else // D++ p
871 sigma *= 4./22.;
872
873 return sigma;
874 }
875 G4double CrossSectionsStrangeness::NDeltaToDeltaSK(Particle const * const p1, Particle const * const p2) {
876 // Nucleon-Delta producing Delta Sigma Kaon cross section
877 //
878 // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
879 //
880 // D++ p (9)
881 // D++ n (15)
882 // D+ p (11)
883 // D+ n (13)
884 //
885 // ratio
886 // D++ p -> S+ K+ D+ (a) (2)
887 // D++ p -> S0 K+ D++ (b) (1)
888 // D++ p -> S+ K0 D++ (c) (6)
889 //
890 // D++ n -> S+ K+ D0 *(d)* (2)
891 // D++ n -> S0 K+ D+ (e) (4)
892 // D++ n -> S- K+ D++ (f) (6)(c)*
893 // D++ n -> S+ K0 D+ (a)* (2)
894 // D++ n -> S0 K0 D++ (b)* (1)*
895 //
896 // D+ p -> S+ K+ D0 (i) (2)*
897 // D+ p -> S0 K+ D+ (j) (1)
898 // D+ p -> S- K+ D++ (k) (2)(g=a)*
899 // D+ p -> S+ K0 D+ (l) (2)
900 // D+ p -> S0 K0 D++ (m) (4)(e)*
901 //
902 // D+ n -> S+ K+ D- *(d)* (2)
903 // D+ n -> S0 K+ D0 (o) (4)
904 // D+ n -> S- K+ D+ (p) (2)*
905 // D+ n -> S+ K0 D0 (i)* (2)*
906 // D+ n -> S0 K0 D+ (j)* (1)*
907 // D+ n -> S- K0 D++ (k)* (2)*
908
909 G4double a = 8.407;
910 G4double b = 2.743;
911 G4double c = 21.18;
912 G4double n_channel = 19.;
913
914 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
915
916 const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
917 const G4double s0 = 8.568E6; // Mev^2
918 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
919
920 if(s <= s0)
921 return 0.;
922
923 G4double sigma = n_channel*a*std::pow(s/s0-1,b)*std::pow(s0/s,c);
924
925 if(iso == 0)// D+ n
926 sigma *= 13./48.;
927 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
928 sigma *= 11./48.;
929 else if (std::abs(iso) == 2)// D++ n
930 sigma *= 15./48.;
931 else // D++ p
932 sigma *= 9./48.;
933
934 return sigma;
935 }
936
937 G4double CrossSectionsStrangeness::NDeltaToNNKKb(Particle const * const p1, Particle const * const p2) {
938 // Nucleon-Delta producing Nucleon-Nucleon Kaon antiKaon cross section
939 //
940 // Total = sigma(NN->NNKKb)*10
941 //
942 // D++ p (6)
943 // D++ n (9)
944 // D+ p (7)
945 // D+ n (8)
946 //
947 // ratio
948 // D++ p -> p p K+ K0b (6)
949 //
950 // D++ n -> p p K+ K- (3)
951 // D++ n -> p p K0 K0b (3)
952 // D++ n -> p n K+ K0b (3)
953 //
954 // D+ p -> p p K+ K- (3)
955 // D+ p -> p p K0 K0b (1)
956 // D+ p -> p n K+ K0b (3)
957 //
958 // D+ n -> p p K0 K- (2)
959 // D+ n -> p n K+ K- (1)
960 // D+ n -> p n K0 K0b (3)
961 // D+ n -> n n K+ K0b (2)
962 //
963
964 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
965
966 G4double sigma = 0.;
967 const G4int iso = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
968 const G4double ener = 0.001*KinematicsUtils::totalEnergyInCM(p1, p2); // GeV
969
970 if(ener <= 2.872)
971 return 0.;
972
973 if(iso == 0)// D+ n
974 sigma = 8* 22./60. * 3. *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
975 else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
976 sigma = 7* 22./60. * 3. *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
977 else if (std::abs(iso) == 2)// D++ n
978 sigma = 9* 22./60. * 3. *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
979 else // D++ p
980 sigma = 6* 22./60. * 3. *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
981
982 return sigma;
983 }
984
985 /// \brief piN to strange cross sections
986
987 G4double CrossSectionsStrangeness::NpiToLK(Particle const * const p1, Particle const * const p2) {
988 //
989 // Pion-Nucleon producing Lambda-Kaon cross sections
990 //
991 // ratio
992 // p pi0 -> L K+ (1/2)
993 // p pi- -> L K0 (1)
994
995 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
996
997 const Particle *pion;
998 const Particle *nucleon;
999 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1000 if(iso == 3 || iso == -3)
1001 return 0.;
1002
1003 if(p1->isPion()){
1004 pion = p1;
1005 nucleon = p2;
1006 }
1007 else{
1008 nucleon = p1;
1009 pion = p2;
1010 }
1011 G4double sigma = 0.;
1012
1013 if(pion->getType() == PiZero)
1014 sigma = 0.5 * p_pimToLK0(pion,nucleon);
1015 else
1016 sigma = p_pimToLK0(pion,nucleon);
1017 return sigma;
1018 }
1019
1020 G4double CrossSectionsStrangeness::p_pimToLK0(Particle const * const p1, Particle const * const p2) {
1021
1022 G4double sigma = 0.;
1023 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1024
1025 if(pLab < 0.911)
1026 return 0.;
1027
1028 sigma = 0.3936*std::pow(pLab,-1.357)-6.052*std::exp(-std::pow(pLab-0.7154,2)/0.02026)-0.16*std::exp(-std::pow(pLab-0.9684,2)/0.001432)+0.489*std::exp(-std::pow(pLab-0.8886,2)/0.08378);
1029 if(sigma < 0.) return 0;
1030 return sigma;
1031 }
1032
1033 G4double CrossSectionsStrangeness::NpiToSK(Particle const * const p1, Particle const * const p2) {
1034 //
1035 // Pion-Nucleon producing Sigma-Kaon cross sections
1036 //
1037 // ratio
1038 // p pi+ (5/3) p pi0 (11/6) p pi- (2)
1039 //
1040 // p pi+ -> S+ K+ (10)
1041 // p pi0 -> S+ K0 (6)*
1042 // p pi0 -> S0 K+ (5)
1043 // p pi- -> S0 K0 (6)*
1044 // p pi- -> S- K+ (6)
1045
1046 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1047
1048 const Particle *pion;
1049 const Particle *nucleon;
1050 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1051
1052 if(p1->isPion()){
1053 pion = p1;
1054 nucleon = p2;
1055 }
1056 else{
1057 nucleon = p1;
1058 pion = p2;
1059 }
1060 G4double sigma = 0.;
1061
1062 if(iso == 3 || iso == -3)
1063 sigma = p_pipToSpKp(pion,nucleon);
1064 else if(pion->getType() == PiZero)
1065 sigma = p_pizToSzKp(pion,nucleon)+p_pimToSzKz(pion,nucleon);
1066 else if(iso == 1 || iso == -1)
1067 sigma = p_pimToSzKz(pion,nucleon)+p_pimToSmKp(pion,nucleon);
1068 else // should never append
1069 sigma = 0.;
1070
1071 return sigma;
1072 }
1073 G4double CrossSectionsStrangeness::p_pimToSmKp(Particle const * const p1, Particle const * const p2) {
1074
1075 G4double sigma = 0.;
1076 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1077
1078 if(pLab < 1.0356)
1079 return 0.;
1080
1081 sigma = 4.352*std::pow(pLab-1.0356,1.006)/(std::pow(pLab+1.0356,0.0978)*std::pow(pLab,5.375));
1082
1083 if(sigma < 0.) // should never append
1084 return 0;
1085
1086 return sigma;
1087 }
1088 G4double CrossSectionsStrangeness::p_pipToSpKp(Particle const * const p1, Particle const * const p2) {
1089
1090 G4double sigma = 0.;
1091 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1092
1093 if(pLab < 1.0428)
1094 return 0.;
1095
1096 sigma = 0.001897*std::pow(pLab-1.0428,2.869)/(std::pow(pLab+1.0428,-16.68)*std::pow(pLab,19.1));
1097
1098 if(sigma < 0.) // should never append
1099 return 0;
1100
1101 return sigma;
1102 }
1103 G4double CrossSectionsStrangeness::p_pizToSzKp(Particle const * const p1, Particle const * const p2) {
1104
1105 G4double sigma = 0.;
1106 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1107
1108 if(pLab < 1.0356)
1109 return 0.;
1110
1111 sigma = 3.624*std::pow(pLab-1.0356,1.4)/std::pow(pLab,5.14);
1112
1113 if(sigma < 0.) // should never append
1114 return 0;
1115
1116 return sigma;
1117 }
1118 G4double CrossSectionsStrangeness::p_pimToSzKz(Particle const * const p1, Particle const * const p2) {
1119
1120 G4double sigma = 0.;
1121 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1122
1123 if((p1->getType() == PiZero && pLab < 1.0356) || (pLab < 1.034))
1124 return 0.;
1125
1126 sigma = 0.3474*std::pow(pLab-1.034,0.07678)/std::pow(pLab,1.627);
1127
1128 if(sigma < 0.) // should never append
1129 return 0;
1130
1131 return sigma;
1132 }
1133
1134 G4double CrossSectionsStrangeness::NpiToLKpi(Particle const * const p1, Particle const * const p2) {
1135 //
1136 // Pion-Nucleon producing Lambda-Kaon-pion cross sections
1137 //
1138 // ratio
1139 // p pi+ (1) p pi0 (3/2) p pi- (2)
1140 //
1141 // p pi0 -> L K+ pi0 (1/2)
1142 // all the others (1)
1143 //
1144 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1145
1146 G4double sigma=0.;
1147 const Particle *pion;
1148 const Particle *nucleon;
1149
1150 if(p1->isPion()){
1151 pion = p1;
1152 nucleon = p2;
1153 }
1154 else{
1155 nucleon = p1;
1156 pion = p2;
1157 }
1158 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1159 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1160
1161 if(pLab < 1.147)
1162 return 0.;
1163
1164 if(iso == 3 || iso == -3)
1165 sigma = 146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1166 else if(pion->getType() == PiZero)
1167 sigma = 1.5*146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1168 else
1169 sigma = 2*146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1170
1171 return sigma;
1172 }
1173
1174 G4double CrossSectionsStrangeness::NpiToSKpi(Particle const * const p1, Particle const * const p2) {
1175 //
1176 // Pion-Nucleon producing Sigma-Kaon-pion cross sections
1177 //
1178 //ratio
1179 // p pi+ (2.25) p pi0 (2.625) p pi-(3)
1180 //
1181 // p pi+ -> S+ pi+ K0 (5/4)
1182 // p pi+ -> S+ pi0 K+ (3/4)
1183 // p pi+ -> S0 pi+ K+ (1/4)
1184 // p pi0 -> S+ pi0 K0 (1/2)
1185 // p pi0 -> S+ pi- K+ (1/2)
1186 // p pi0 -> S0 pi+ K0 (3/4)
1187 // p pi0 -> S0 pi0 K+ (3/8)
1188 // p pi0 -> S- pi+ K+ (1/2)
1189 // p pi- -> S+ pi- K0 (3/8)
1190 // p pi- -> S0 pi0 K0 (5/8)
1191 // p pi- -> S0 pi- K+ (5/8)
1192 // p pi- -> S- pi+ K0 (1)
1193 // p pi- -> S- pi0 K+ (3/8)
1194
1195 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1196
1197 G4double sigma=0.;
1198 const Particle *pion;
1199 const Particle *nucleon;
1200
1201 if(p1->isPion()){
1202 pion = p1;
1203 nucleon = p2;
1204 }
1205 else{
1206 nucleon = p1;
1207 pion = p2;
1208 }
1209 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1210 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1211
1212 if(pLab <= 1.3041)
1213 return 0.;
1214
1215 if(iso == 3 || iso == -3)
1216 sigma = 2.25*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1217 else if(pion->getType() == PiZero)
1218 sigma = 2.625*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1219 else
1220 sigma = 3.*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1221
1222 return sigma;
1223 }
1224
1225 G4double CrossSectionsStrangeness::NpiToLK2pi(Particle const * const p1, Particle const * const p2) {
1226 //
1227 // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1228 //
1229 // p pi+ (2) p pi0 (1.75) p pi- (2.5)
1230 //
1231 // p pi+ -> L K+ pi+ pi0 (1)
1232 // p pi+ -> L K0 pi+ pi+ (1)
1233 // p pi0 -> L K+ pi0 pi0 (1/4)
1234 // p pi0 -> L K+ pi+ pi- (1)
1235 // p pi0 -> L K0 pi+ pi0 (1/2)
1236 // p pi- -> L K+ pi0 pi- (1)
1237 // p pi- -> L K0 pi+ pi- (1)
1238 // p pi- -> L K0 pi0 pi0 (1/2)
1239
1240 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1241
1242 G4double sigma=0.;
1243 const Particle *pion;
1244 const Particle *nucleon;
1245
1246 if(p1->isPion()){
1247 pion = p1;
1248 nucleon = p2;
1249 }
1250 else{
1251 nucleon = p1;
1252 pion = p2;
1253 }
1254 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1255 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1256
1257 if(pLab <= 1.4162)
1258 return 0.;
1259
1260 if(iso == 3 || iso == -3)
1261 sigma = 2*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1262 else if(pion->getType() == PiZero)
1263 sigma = 1.75*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1264 else
1265 sigma = 2.5*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1266
1267 return sigma;
1268 }
1269
1270 G4double CrossSectionsStrangeness::NpiToSK2pi(Particle const * const p1, Particle const * const p2) {
1271 //
1272 // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1273 //
1274 // ratio
1275 // p pi+ (3.25) p pi0 (3.5) p pi- (3.75)
1276 //
1277 // p pi+ -> S+ K+ pi+ pi- (1)
1278 // p pi+ -> S+ K+ pi0 pi0 (1/4)
1279 // p pi+ -> S0 K+ pi+ pi0 (1/2)
1280 // p pi+ -> S- K+ pi+ pi+ (1/4)
1281 // p pi+ -> S+ K0 pi+ pi0 (1)
1282 // p pi+ -> S0 K0 pi+ pi+ (1/4)
1283 //
1284 // p pi0 -> S+ K+ pi0 pi- (1/2)
1285 // p pi0 -> S0 K+ pi+ pi- (1/2)
1286 // p pi0 -> S0 K+ pi0 pi0 (1/4)
1287 // p pi0 -> S- K+ pi+ pi0 (1/4)
1288 // p pi0 -> S+ K0 pi+ pi- (1)
1289 // p pi0 -> S+ K0 pi0 pi0 (1/4)
1290 // p pi0 -> S0 K0 pi+ pi0 (1/4)
1291 // p pi0 -> S- K0 pi+ pi+ (1/2)
1292 //
1293 // p pi- -> S+ K+ pi- pi- (1/4)
1294 // p pi- -> S0 K+ pi0 pi- (1/2)
1295 // p pi- -> S- K+ pi+ pi- (1/4)
1296 // p pi- -> S- K+ pi0 pi0 (1/4)
1297 // p pi- -> S+ K0 pi0 pi- (1/2)
1298 // p pi- -> S0 K0 pi+ pi- (1)
1299 // p pi- -> S0 K0 pi0 pi0 (1/2)
1300 // p pi- -> S- K0 pi+ pi0 (1/2)
1301
1302 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1303
1304 G4double sigma=0.;
1305 const Particle *pion;
1306 const Particle *nucleon;
1307
1308 if(p1->isPion()){
1309 pion = p1;
1310 nucleon = p2;
1311 }
1312 else{
1313 nucleon = p1;
1314 pion = p2;
1315 }
1316 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1317 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1318
1319 if(pLab <= 1.5851)
1320 return 0.;
1321
1322 if(iso == 3 || iso == -3)
1323 sigma = 3.25*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1324 else if(pion->getType() == PiZero)
1325 sigma = 3.5*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1326 else
1327 sigma = 3.75*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1328
1329 return sigma;
1330 }
1331
1332 G4double CrossSectionsStrangeness::NpiToNKKb(Particle const * const p1, Particle const * const p2) {
1333 //
1334 // Pion-Nucleon producing Nucleon-Kaon-antiKaon cross sections
1335 //
1336 // ratio
1337 // p pi+ (1/2) p pi0 (3/2) p pi- (5/2)
1338 //
1339 // p pi+ -> p K+ K0b (1/2)
1340 // p pi0 -> p K0 K0b (1/4)
1341 // p pi0 -> p K+ K- (1/4)
1342 // p pi0 -> n K+ K0b (1)
1343 // p pi- -> p K0 K- (1/2)
1344 // p pi- -> n K+ K- (1)
1345 // p pi- -> n K0 K0b (1)
1346
1347 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1348
1349 const Particle *particle1;
1350 const Particle *particle2;
1351
1352 if(p1->isPion()){
1353 particle1 = p1;
1354 particle2 = p2;
1355 }
1356 else{
1357 particle1 = p2;
1358 particle2 = p1;
1359 }
1360
1361 G4double sigma = 0.;
1362
1363 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1364
1365 if(particle1->getType() == PiZero){
1366 if(pLab < 1.5066) return 0.;
1367 else if(pLab < 30.) sigma = 3./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1368 else return 0.;
1369 }
1370 else if((particle1->getType() == PiPlus && particle2->getType() == Neutron) || (particle1->getType() == PiMinus && particle2->getType() == Proton)){
1371 if(pLab < 1.5066) return 0.;
1372 else if(pLab < 30.) sigma = 5./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1373 else return 0.;
1374 }
1375 else{
1376 if(pLab < 1.5066) return 0.;
1377 else if(pLab < 30.) sigma = 1./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1378 else return 0.;
1379 }
1380 return sigma;
1381 }
1382
1383 G4double CrossSectionsStrangeness::NpiToMissingStrangeness(Particle const * const p1, Particle const * const p2) {
1384 //
1385 // Pion-Nucleon missing strangeness production cross sections
1386 //
1387 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1388
1389 const Particle *pion;
1390 const Particle *nucleon;
1391
1392 if(p1->isPion()){
1393 pion = p1;
1394 nucleon = p2;
1395 }
1396 else{
1397 pion = p2;
1398 nucleon = p1;
1399 }
1400
1401 G4double sigma = 0.;
1402
1403 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(pion,nucleon); // GeV
1404 if(pLab < 2.2) return 0.;
1405
1406 if(pion->getType() == PiZero){
1407 if(pLab < 30.) sigma = 4.4755*std::pow((pLab - 2.2),1.927)/std::pow(pLab,1.89343);
1408 else return 0.;
1409 }
1410 else if((pion->getType() == PiPlus && nucleon->getType() == Neutron) || (pion->getType() == PiMinus && nucleon->getType() == Proton)){
1411 if(pLab < 30.) sigma = 5.1*std::pow((pLab - 2.2),1.854)/std::pow(pLab,1.904);
1412 else return 0.;
1413 }
1414 else{
1415 if(pLab < 30.) sigma = 3.851*std::pow((pLab - 2.2),2)/std::pow(pLab,1.88286);
1416 else return 0.;
1417 }
1418 return sigma;
1419 }
1420
1421 /// \brief NY cross sections
1422
1423 G4double CrossSectionsStrangeness::NLToNS(Particle const * const p1, Particle const * const p2) {
1424 //
1425 // Nucleon-Lambda producing Nucleon-Sigma cross sections
1426 //
1427 // ratio
1428 // p L -> p S0 (1/2)
1429 // p L -> n S+ (1)
1430
1431
1432 // assert((p1->isLambda() && p2->isNucleon()) || (p2->isLambda() && p1->isNucleon()));
1433
1434 G4double sigma = 0.;
1435
1436 const Particle *particle1;
1437 const Particle *particle2;
1438
1439 if(p1->isLambda()){
1440 particle1 = p1;
1441 particle2 = p2;
1442 }
1443 else{
1444 particle1 = p2;
1445 particle2 = p1;
1446 }
1447
1448 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2);
1449
1450 if(pLab < 0.664)
1451 return 0.;
1452
1453 sigma = 3 * 8.74*std::pow((pLab-0.664),0.438)/std::pow(pLab,2.717); // 3 * L p -> S0 p
1454
1455 return sigma;
1456 }
1457
1458 G4double CrossSectionsStrangeness::NSToNL(Particle const * const p1, Particle const * const p2) {
1459 //
1460 // Nucleon-Lambda producing Nucleon-Sigma cross sections
1461 //
1462 // ratio
1463 // p S0 -> p L (1/2)
1464 // p S- -> n L (1)
1465
1466 // assert((p1->isSigma() && p2->isNucleon()) || (p2->isSigma() && p1->isNucleon()));
1467
1468 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1469 if(iso == 3 || iso == -3)
1470 return 0.;
1471
1472 G4double sigma;
1473 const Particle *particle1;
1474 const Particle *particle2;
1475
1476 if(p1->isSigma()){
1477 particle1 = p1;
1478 particle2 = p2;
1479 }
1480 else{
1481 particle1 = p2;
1482 particle2 = p1;
1483 }
1484 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1485
1486 if(particle1->getType() == SigmaZero){
1487 if(pLab < 0.1) return 100.; // cut-max
1488 sigma = 8.23*std::pow(pLab,-1.087);
1489 }
1490 else{
1491 if(pLab < 0.1) return 200.; // cut-max
1492 sigma = 16.46*std::pow(pLab,-1.087);
1493 }
1494 return sigma;
1495 }
1496
1497 G4double CrossSectionsStrangeness::NSToNS(Particle const * const p1, Particle const * const p2) {
1498
1499 // assert((p1->isSigma() && p2->isNucleon()) || (p2->isSigma() && p1->isNucleon()));
1500
1501 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1502 if(iso == 3 || iso == -3)
1503 return 0.; // only quasi-elastic here
1504
1505 const Particle *particle1;
1506 const Particle *particle2;
1507
1508 if(p1->isSigma()){
1509 particle1 = p1;
1510 particle2 = p2;
1511 }
1512 else{
1513 particle1 = p2;
1514 particle2 = p1;
1515 }
1516 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1517
1518 if(particle2->getType() == Neutron && pLab < 0.162) return 0.;
1519 else if(pLab < 0.1035) return 200.; // cut-max
1520
1521 return 13.79*std::pow(pLab,-1.181);
1522 }
1523
1524 /// \brief NK cross sections
1525
1526 G4double CrossSectionsStrangeness::NKToNK(Particle const * const p1, Particle const * const p2) {
1527 //
1528 // Nucleon-Kaon quasi-elastic cross sections
1529 //
1530 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1531
1532 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1533
1534 if(iso != 0)
1535 return 0.;
1536
1537 const Particle *particle1;
1538 const Particle *particle2;
1539
1540 if(p1->isKaon()){
1541 particle1 = p1;
1542 particle2 = p2;
1543 }
1544 else{
1545 particle1 = p2;
1546 particle2 = p1;
1547 }
1548
1549 G4double sigma = 0.;
1550 G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1551
1552 if(particle1->getType() == Proton)
1553 pLab += 2*0.0774;
1554
1555 if(pLab <= 0.0774)
1556 return 0.;
1557
1558 sigma = 12.84*std::pow((pLab-0.0774),18.19)/std::pow((pLab),20.41);
1559
1560 return sigma;
1561 }
1562
1563 G4double CrossSectionsStrangeness::NKToNKpi(Particle const * const p1, Particle const * const p2) {
1564 //
1565 // Nucleon-Kaon producing Nucleon-Kaon-pion cross sections
1566 //
1567 // Ratio determined by meson symmetry using only "resonante" diagram (with Delta or K*)
1568 //
1569 // ratio: K+ p (5) K0 p (5.545)
1570 //
1571 // K+ p -> p K+ pi0 1.2
1572 // K+ p -> p K0 pi+ 3
1573 // K+ p -> n K+ pi+ 0.8
1574 // K0 p -> p K+ pi- 1
1575 // K0 p -> p K0 pi0 0.845
1576 // K0 p -> n K+ pi0 1.47
1577 // K0 p -> n K0 pi+ 2.23
1578 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1579
1580 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1581
1582 const Particle *particle1;
1583 const Particle *particle2;
1584
1585 if(p1->isKaon()){
1586 particle1 = p1;
1587 particle2 = p2;
1588 }
1589 else{
1590 particle1 = p2;
1591 particle2 = p1;
1592 }
1593
1594 G4double sigma = 0.;
1595 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1596
1597 if(pLab <= 0.53)
1598 return 0.;
1599
1600 if(iso == 0)
1601 sigma = 5.55*116.8*std::pow(pLab-0.53,6.874)/std::pow(pLab,10.11);
1602 else
1603 sigma = 5.*116.8*std::pow(pLab-0.53,6.874)/std::pow(pLab,10.11);;
1604 return sigma;
1605 }
1606
1607 G4double CrossSectionsStrangeness::NKToNK2pi(Particle const * const p1, Particle const * const p2) {
1608 //
1609 // Nucleon-Kaon producing Nucleon-Kaon-2pion cross sections
1610 //
1611 // p K+ (2.875) p K0 (3.125)
1612 //
1613 // p K+ -> p K+ pi+ pi- (1)
1614 // p K+ -> p K+ pi0 pi0 (1/8)
1615 // p K+ -> p K0 pi+ pi0 (1)
1616 // p K+ -> n K+ pi+ pi0 (1/2)
1617 // p K+ -> n K0 pi+ pi+ (1/4)
1618 // p K0 -> p K+ pi0 pi- (1)
1619 // p K0 -> p K0 pi+ pi- (1)
1620 // p K0 -> p K0 pi0 pi0 (1/8)
1621 // p K0 -> n K+ pi+ pi- (1/4)
1622 // p K0 -> n K+ pi0 pi0 (1/4)
1623 // p K0 -> n K0 pi+ pi0 (1/2)
1624
1625 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1626
1627 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1628
1629 const Particle *particle1;
1630 const Particle *particle2;
1631
1632 if(p1->isKaon()){
1633 particle1 = p1;
1634 particle2 = p2;
1635 }
1636 else{
1637 particle1 = p2;
1638 particle2 = p1;
1639 }
1640
1641 G4double sigma = 0.;
1642 const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1643
1644 if(pLab < 0.812)
1645 sigma = 0.;
1646 else if(pLab < 1.744)
1647 sigma = 26.41*std::pow(pLab-0.812,7.138)/std::pow(pLab,5.337);
1648 else if(pLab < 3.728)
1649 sigma = 1572.*std::pow(pLab-0.812,9.069)/std::pow(pLab,12.44);
1650 else
1651 sigma = 60.23*std::pow(pLab-0.812,5.084)/std::pow(pLab,6.72);
1652
1653 if(iso == 0)
1654 sigma *= 3.125;
1655 else
1656 sigma *= 2.875;
1657
1658 return sigma;
1659 }
1660
1661 /// \brief NKb cross sections
1662
1663 G4double CrossSectionsStrangeness::NKbToNKb(Particle const * const p1, Particle const * const p2) {
1664 //
1665 // Nucleon-antiKaon quasi-elastic cross sections
1666 //
1667 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1668
1669 G4double sigma=0.;
1670 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1671
1672 const Particle *antikaon;
1673 const Particle *nucleon;
1674
1675 if (p1->isAntiKaon()) {
1676 antikaon = p1;
1677 nucleon = p2;
1678 }
1679 else {
1680 antikaon = p2;
1681 nucleon = p1;
1682 }
1683
1684 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1685
1686 if(iso != 0) // K0b p and K- n -> forbidden: quasi-elastic diffusion only
1687 return 0;
1688 else if(nucleon->getType() == Proton){ // K- p -> K0b n
1689 if(pLab < 0.08921)
1690 return 0.;
1691 else if(pLab < 0.2)
1692 sigma = 0.4977*std::pow(pLab - 0.08921,0.5581)/std::pow(pLab,2.704);
1693 else if(pLab < 0.73)
1694 sigma = 2.*std::pow(pLab,-1.2) + 6.493*std::exp(-0.5*std::pow((pLab-0.3962)/0.02,2));
1695 else if(pLab < 1.38)
1696 sigma = 2.3*std::pow(pLab,-0.9) + 1.1*std::exp(-0.5*std::pow((pLab-0.82)/0.04,2)) + 5.*std::exp(-0.5*std::pow((pLab-1.04)/0.1,2));
1697 else if(pLab < 30.)
1698 sigma = 2.5*std::pow(pLab,-1.68) + 0.7*std::exp(-0.5*std::pow((pLab-1.6)/0.2,2)) + 0.2*std::exp(-0.5*std::pow((pLab-2.3)/0.2,2));
1699 else sigma = 0.;
1700 }
1701 else{ // K0b n -> K- p (same as K- p but without threshold)
1702 if(pLab < 0.1)
1703 sigma = 30.;
1704 else if(pLab < 0.73)
1705 sigma = 2.*std::pow(pLab,-1.2) + 6.493*std::exp(-0.5*std::pow((pLab-0.3962)/0.02,2));
1706 else if(pLab < 1.38)
1707 sigma = 2.3*std::pow(pLab,-0.9) + 1.1*std::exp(-0.5*std::pow((pLab-0.82)/0.04,2)) + 5.*std::exp(-0.5*std::pow((pLab-1.04)/0.1,2));
1708 else if(pLab < 30.)
1709 sigma = 2.5*std::pow(pLab,-1.68) + 0.7*std::exp(-0.5*std::pow((pLab-1.6)/0.2,2)) + 0.2*std::exp(-0.5*std::pow((pLab-2.3)/0.2,2));
1710 else sigma = 0.;
1711 }
1712 return sigma;
1713 }
1714
1715 G4double CrossSectionsStrangeness::NKbToSpi(Particle const * const p1, Particle const * const p2) {
1716 //
1717 // Nucleon-antiKaon producing Sigma-pion cross sections
1718 //
1719 // ratio
1720 // p K0b (4/3) p K- (13/6)
1721 //
1722 // p K0b -> S+ pi0 (2/3)
1723 // p K0b -> S0 pi+ (2/3)
1724 // p K- -> S+ pi- (1)
1725 // p K- -> S0 pi0 (1/2)
1726 // p K- -> S- pi+ (2/3)
1727
1728 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1729
1730 G4double sigma=0.;
1731 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1732
1733 const Particle *antikaon;
1734 const Particle *nucleon;
1735
1736 if (p1->isAntiKaon()) {
1737 antikaon = p1;
1738 nucleon = p2;
1739 }
1740 else {
1741 antikaon = p2;
1742 nucleon = p1;
1743 }
1744
1745 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1746
1747 if(iso == 0){
1748 if(pLab < 0.1)
1749 return 152.0; // 70.166*13/6
1750 else
1751 sigma = 13./6.*(1.4*std::pow(pLab,-1.7)+1.88*std::exp(-std::pow(pLab-0.747,2)/0.005)+8*std::exp(-std::pow(pLab-0.4,2)/0.002)+0.8*std::exp(-std::pow(pLab-1.07,2)/0.01));
1752 }
1753 else{
1754 if(pLab < 0.1)
1755 return 93.555; // 70.166*4/3
1756 else
1757 sigma = 4./3.*(1.4*std::pow(pLab,-1.7)+1.88*std::exp(-std::pow(pLab-0.747,2)/0.005)+8*std::exp(-std::pow(pLab-0.4,2)/0.002)+0.8*std::exp(-std::pow(pLab-1.07,2)/0.01));
1758 //sigma = 4./3.*(1.4*std::pow(pLab,-1.7)+1.88*std::exp(-std::pow(pLab-0.747,2)/0.005)+0.8*std::exp(-std::pow(pLab-1.07,2)/0.01));
1759 }
1760
1761 return sigma;
1762 }
1763
1764 G4double CrossSectionsStrangeness::NKbToLpi(Particle const * const p1, Particle const * const p2) {
1765 //
1766 // Nucleon-antiKaon producing Lambda-pion cross sections
1767 //
1768 // ratio
1769 // p K0b (1) p K- (1/2)
1770 //
1771 // p K- -> L pi0 (1/2)
1772 // p K0b -> L pi+ (1)
1773 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1774
1775 G4double sigma = 0.;
1776 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1777
1778 const Particle *antikaon;
1779 const Particle *nucleon;
1780
1781 if (p1->isAntiKaon()) {
1782 antikaon = p1;
1783 nucleon = p2;
1784 }
1785 else {
1786 antikaon = p2;
1787 nucleon = p1;
1788 }
1789 if(iso == 0)
1790 sigma = p_kmToL_pz(antikaon,nucleon);
1791 else
1792 sigma = 2*p_kmToL_pz(antikaon,nucleon);
1793
1794 return sigma;
1795 }
1796 G4double CrossSectionsStrangeness::p_kmToL_pz(Particle const * const p1, Particle const * const p2) {
1797 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(p1, p2); // GeV
1798 G4double sigma = 0.;
1799 if(pLab < 0.086636)
1800 sigma = 40.24;
1801 else if(pLab < 0.5)
1802 sigma = 0.97*std::pow(pLab,-1.523);
1803 else if(pLab < 2.)
1804 sigma = 1.23*std::pow(pLab,-1.467)+0.872*std::exp(-std::pow(pLab-0.749,2)/0.0045)+2.337*std::exp(-std::pow(pLab-0.957,2)/0.017)+0.476*std::exp(-std::pow(pLab-1.434,2)/0.136);
1805 else if(pLab < 30.)
1806 sigma = 3.*std::pow(pLab,-2.57);
1807 else
1808 sigma = 0.;
1809
1810 return sigma;
1811 }
1812
1813 G4double CrossSectionsStrangeness::NKbToS2pi(Particle const * const p1, Particle const * const p2) {
1814 //
1815 // Nucleon-antiKaon producing Sigma-2pion cross sections
1816 //
1817 // ratio
1818 // p K0b (29/12) p K- (59/24)
1819 //
1820 // p K0b -> S+ pi+ pi- (2/3)
1821 // p K0b -> S+ pi0 pi0 (1/4)
1822 // p K0b -> S0 pi+ pi0 (5/6)
1823 // p K0b -> S- pi+ pi+ (2/3)
1824 // p K- -> S+ pi0 pi- (1)
1825 // p K- -> S0 pi+ pi- (2/3)
1826 // p K- -> S0 pi0 pi0 (1/8)
1827 // p K- -> S- pi+ pi0 (2/3)
1828
1829 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1830
1831 G4double sigma=0.;
1832 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1833
1834 const Particle *antikaon;
1835 const Particle *nucleon;
1836
1837 if (p1->isAntiKaon()) {
1838 antikaon = p1;
1839 nucleon = p2;
1840 }
1841 else {
1842 antikaon = p2;
1843 nucleon = p1;
1844 }
1845
1846 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1847
1848 if(pLab < 0.260)
1849 return 0.;
1850
1851 if(iso == 0)
1852 sigma = 29./12.*3./2.*(49.96*std::pow(pLab-0.260,6.398)/std::pow(pLab+0.260,9.732)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));
1853 else
1854 sigma = 54./24.*3./2.*(49.96*std::pow(pLab-0.260,6.398)/std::pow(pLab+0.260,9.732)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));
1855
1856 /*
1857 if(iso == 0)
1858 sigma = 29./12.*(85.46*std::pow(pLab-0.226,8.118)/std::pow(pLab+0.226,11.69)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));
1859 else
1860 sigma = 54./24.*(85.46*std::pow(pLab-0.226,8.118)/std::pow(pLab+0.226,11.69)+0.1451*std::exp(-std::pow(pLab-0.4031,2)/0.00115));*/
1861
1862 return sigma;
1863 }
1864
1865 G4double CrossSectionsStrangeness::NKbToL2pi(Particle const * const p1, Particle const * const p2) {
1866 //
1867 // Nucleon-antiKaon producing Lambda-2pion cross sections
1868 //
1869 // ratio
1870 // p K0b -> L pi+ pi0 (1)
1871 // p K- -> L pi+ pi- (1)
1872 // p K- -> L pi0 pi0 (1/4)
1873
1874 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1875
1876 G4double sigma = 0.;
1877 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1878
1879 const Particle *antikaon;
1880 const Particle *nucleon;
1881
1882 if (p1->isAntiKaon()) {
1883 antikaon = p1;
1884 nucleon = p2;
1885 }
1886 else {
1887 antikaon = p2;
1888 nucleon = p1;
1889 }
1890
1891 if(iso == 0)
1892 sigma = 1.25*p_kmToL_pp_pm(antikaon,nucleon);
1893 else
1894 sigma = p_kmToL_pp_pm(antikaon,nucleon);
1895
1896 return sigma;
1897 }
1898 G4double CrossSectionsStrangeness::p_kmToL_pp_pm(Particle const * const p1, Particle const * const p2) {
1899 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(p1, p2); // GeV
1900 G4double sigma = 0.;
1901 if(pLab < 0.97)
1902 sigma = 6364.*std::pow(pLab,6.07)/std::pow(pLab+1.,10.58)+2.158*std::exp(-std::pow((pLab-0.395)/.01984,2)/2.);
1903 else if(pLab < 30)
1904 sigma = 46.3*std::pow(pLab,0.62)/std::pow(pLab+1.,3.565);
1905 else
1906 sigma = 0.;
1907
1908 return sigma;
1909 }
1910
1911 G4double CrossSectionsStrangeness::NKbToNKbpi(Particle const * const p1, Particle const * const p2) {
1912 //
1913 // Nucleon-antiKaon producing Nucleon-antiKaon-pion cross sections
1914 //
1915 // ratio
1916 // p K- (28) p K0b (20)
1917 //
1918 // p K- -> p K- pi0 (6)*
1919 // p K- -> p K0b pi- (7)*
1920 // p K- -> n K- pi+ (9)*
1921 // p K- -> n K0b pi0 (6)
1922 // p K0b -> p K0b pi0 (4)
1923 // p K0b -> p K- pi+ (10)*
1924 // p K0b -> n K0b pi+ (6)*
1925 //
1926
1927 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1928
1929 G4double sigma=0.;
1930 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1931
1932 const Particle *antikaon;
1933 const Particle *nucleon;
1934
1935 if (p1->isAntiKaon()) {
1936 antikaon = p1;
1937 nucleon = p2;
1938 }
1939 else {
1940 antikaon = p2;
1941 nucleon = p1;
1942 }
1943
1944 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1945
1946 if(pLab < 0.526)
1947 return 0.;
1948
1949 if(iso == 0)
1950 sigma = 28. * 10.13*std::pow(pLab-0.526,5.846)/std::pow(pLab,8.343);
1951 else
1952 sigma = 20. * 10.13*std::pow(pLab-0.526,5.846)/std::pow(pLab,8.343);
1953
1954 return sigma;
1955 }
1956
1957 G4double CrossSectionsStrangeness::NKbToNKb2pi(Particle const * const p1, Particle const * const p2) {
1958 //
1959 // Nucleon-antiKaon producing Nucleon-antiKaon-2pion cross sections
1960 //
1961 // ratio
1962 // p K0b (4.25) p K- (4.75)
1963 //
1964 // p K0b -> p K0b pi+ pi- (1)
1965 // p K0b -> p K0b pi0 pi0 (1/4)
1966 // p K0b -> p K- pi+ pi0 (1)
1967 // p K0b -> n K0b pi+ pi0 (1)
1968 // p K0b -> n K- pi+ pi+ (1)
1969 // p K- -> p K0b pi0 pi- (1)
1970 // p K- -> p K- pi+ pi- (1)
1971 // p K- -> p K- pi0 pi0 (1/4)
1972 // p K- -> n K0b pi+ pi- (1)
1973 // p K- -> n K0b pi0 pi0 (1/2)
1974 // p K- -> n K- pi+ pi0 (1)
1975
1976 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1977
1978 G4double sigma=0.;
1979 const G4int iso=ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1980
1981 const Particle *antikaon;
1982 const Particle *nucleon;
1983
1984 if (p1->isAntiKaon()) {
1985 antikaon = p1;
1986 nucleon = p2;
1987 }
1988 else {
1989 antikaon = p2;
1990 nucleon = p1;
1991 }
1992
1993 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1994
1995 if(pLab < 0.85)
1996 return 0.;
1997
1998 if(iso == 0)
1999 sigma = 4.75 * 26.8*std::pow(pLab-0.85,4.9)/std::pow(pLab,6.34);
2000 else
2001 sigma = 4.25 * 26.8*std::pow(pLab-0.85,4.9)/std::pow(pLab,6.34);
2002
2003 return sigma;
2004 }
2005
2006
2007 } // namespace G4INCL
2008
2009