ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCrossSectionsStrangeness.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCrossSectionsStrangeness.cc
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 
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 
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 
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 
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 
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 
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 
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 
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 
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