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