Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 
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;
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 
65 
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())) {
95  } else if((p1->isNucleon() && p2->isOmega()) ||
96  (p1->isOmega() && p2->isNucleon())) {
98  } else if((p1->isNucleon() && p2->isEtaPrime()) ||
99  (p1->isEtaPrime() && p2->isNucleon())) {
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())){
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);
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;
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;
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;
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;
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.;
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 
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
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  // D++ n -> p L K+ (3)
738  //
739  // D+ p -> p L K+ (1)
740  //
741  // D+ n -> p L K0 (1)
742  // D+ n -> n L K+ (1)
743  //return 0.;
744 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
745 
747  if(std::abs(iso) == 4) return 0.;
748 
749  G4double sigma = 0.;
750 
751  const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // MeV^2
752  const G4double s0 = 6.511E6; // MeV^2
753 
754  if(s <= s0) return 0.;
755 
756  sigma = 4.*4.169*std::pow(s/s0-1,2.227)*std::pow(s0/s,2.511);
757 
758  //const G4double pLab = sdt::sqrt(s*s/(4*ParticleTable::effectiveNucleonMass2)-s)*0.001;
759  //sigma = 3*1.11875*std::pow((pLab-2.3508),1.0951)/std::pow((pLab+2.3508),2.0958); // NDelta sim to NN
760 
761  if(iso == 0){ // D+ n
762  sigma *= 2./6.;
763  }
764  else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType())){// D+ p
765  sigma *= 1./6.;
766  }
767  else{// D++ n
768  sigma *= 3./6.;
769  }
770  return sigma;
771  }
772  G4double CrossSectionsStrangeness::NDeltaToNSK(Particle const * const p1, Particle const * const p2) {
773  // Nucleon-Delta producing Nucleon Sigma Kaon cross section
774  //
775  // 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)
776  //
777  // D++ p -> p S+ K+ (6)
778  //
779  // D++ n -> p S+ K0 (3) ****
780  // D++ n -> p S0 K+ (3)
781  // D++ n -> n S+ K+ (3)
782  //
783  // D+ p -> p S+ K0 (2)
784  // D+ p -> p S0 K+ (2)
785  // D+ p -> n S+ K+ (3)
786  //
787  // D+ n -> p S0 K0 (3)
788  // D+ n -> p S- K+ (2)
789  // D+ n -> n S+ K0 (2)
790  // D+ n -> n S0 K+ (2)
791 
792 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
793 
794  G4double sigma = 0.;
795 
796  const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
797  const G4double s0 = 6.935E6; // Mev^2
799 
800  if(s <= s0)
801  return 0.;
802 
803  sigma = 11.*39.54*std::pow(s/s0-1,2.799)*std::pow(s0/s,6.303);
804 
805  //const G4double pLab = sdt::sqrt(s*s/(4*ParticleTable::effectiveNucleonMass2)-s)*0.001;
806  //sigma = 22./12./2. * 4.75*6.38*std::pow(pLab-2.593,2.1)/std::pow(pLab,4.162); // NDelta sim to NN
807 
808  if(iso == 0)// D+ n
809  sigma *= 9./31.;
810  else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
811  sigma *= 7./31.;
812  else if (std::abs(iso) == 2)// D++ n
813  sigma *= 9./31.;
814  else // D++ p
815  sigma *= 6./31.;
816 
817  return sigma;
818  }
820  // Nucleon-Delta producing Delta Lambda Kaon cross section
821  //
822  // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
823  //
824  // ratio
825  // D++ p -> L K+ D++ (4)
826  //
827  // D++ n -> L K+ D+ (3)
828  // D++ n -> L K0 D++ (4)
829  //
830  // D+ p -> L K0 D++ (3)
831  // D+ p -> L K+ D+ (2)
832  //
833  // D+ n -> L K+ D0 (4)
834  // D+ n -> L K0 D+ (2)
835  //return 0.;
836 
837 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
838 
839  const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
840  const G4double s0 = 8.096E6; // Mev^2
842 
843  if(s <= s0)
844  return 0.;
845 
846  G4double sigma = 7.*2.679*std::pow(s/s0-1,2.280)*std::pow(s0/s,5.086);
847 
848  if(iso == 0)// D+ n
849  sigma *= 6./22.;
850  else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
851  sigma *= 5./22.;
852  else if (std::abs(iso) == 2)// D++ n
853  sigma *= 7./22.;
854  else // D++ p
855  sigma *= 4./22.;
856 
857  return sigma;
858  }
860  // Nucleon-Delta producing Delta Sigma Kaon cross section
861  //
862  // XS from K. Tsushima, A. Sibirtsev, A. W. Thomas, and G. Q. Li. Phys.Rev.C 59, 369
863  //
864  // D++ p (9)
865  // D++ n (15)
866  // D+ p (11)
867  // D+ n (13)
868  //
869  // ratio
870  // D++ p -> S+ K+ D+ (a) (2)
871  // D++ p -> S0 K+ D++ (b) (1)
872  // D++ p -> S+ K0 D++ (c) (6)
873  //
874  // D++ n -> S+ K+ D0 *(d)* (2)
875  // D++ n -> S0 K+ D+ (e) (4)
876  // D++ n -> S- K+ D++ (f) (6)(c)*
877  // D++ n -> S+ K0 D+ (a)* (2)
878  // D++ n -> S0 K0 D++ (b)* (1)*
879  //
880  // D+ p -> S+ K+ D0 (i) (2)*
881  // D+ p -> S0 K+ D+ (j) (1)
882  // D+ p -> S- K+ D++ (k) (2)(g=a)*
883  // D+ p -> S+ K0 D+ (l) (2)
884  // D+ p -> S0 K0 D++ (m) (4)(e)*
885  //
886  // D+ n -> S+ K+ D- *(d)* (2)
887  // D+ n -> S0 K+ D0 (o) (4)
888  // D+ n -> S- K+ D+ (p) (2)*
889  // D+ n -> S+ K0 D0 (i)* (2)*
890  // D+ n -> S0 K0 D+ (j)* (1)*
891  // D+ n -> S- K0 D++ (k)* (2)*
892  //return 0.;
893 
894 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
895 
896  const G4double s = KinematicsUtils::squareTotalEnergyInCM(p1 ,p2); // Mev^^2
897  const G4double s0 = 8.568E6; // Mev^2
899 
900  if(s <= s0)
901  return 0.;
902 
903  G4double sigma = 19.*21.18*std::pow(s/s0-1,2.743)*std::pow(s0/s,8.407);
904 
905  if(iso == 0)// D+ n
906  sigma *= 13./48.;
907  else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
908  sigma *= 11./48.;
909  else if (std::abs(iso) == 2)// D++ n
910  sigma *= 15./48.;
911  else // D++ p
912  sigma *= 9./48.;
913 
914  return sigma;
915  }
916 
918  // Nucleon-Delta producing Nucleon-Nucleon Kaon antiKaon cross section
919  //
920  // Total = sigma(NN->NNKKb)*10
921  //
922  // D++ p (6)
923  // D++ n (9)
924  // D+ p (7)
925  // D+ n (8)
926  //
927  // ratio
928  // D++ p -> p p K+ K0b (6)
929  //
930  // D++ n -> p p K+ K- (3)
931  // D++ n -> p p K0 K0b (3)
932  // D++ n -> p n K+ K0b (3)
933  //
934  // D+ p -> p p K+ K- (3)
935  // D+ p -> p p K0 K0b (1)
936  // D+ p -> p n K+ K0b (3)
937  //
938  // D+ n -> p p K0 K- (2)
939  // D+ n -> p n K+ K- (1)
940  // D+ n -> p n K0 K0b (3)
941  // D+ n -> n n K+ K0b (2)
942  //
943 
944 // assert((p1->isNucleon() && p2->isResonance()) || (p2->isNucleon() && p1->isResonance()));
945 
946  G4double sigma = 0.;
948  const G4double ener = 0.001*KinematicsUtils::totalEnergyInCM(p1, p2); // GeV
949 
950  if(ener <= 2.872)
951  return 0.;
952 
953  if(iso == 0)// D+ n
954  sigma = 8* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
955  else if (ParticleTable::getIsospin(p1->getType()) == ParticleTable::getIsospin(p2->getType()))// D+ p
956  sigma = 7* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
957  else if (std::abs(iso) == 2)// D++ n
958  sigma = 9* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
959  else // D++ p
960  sigma = 6* 14. * 5./19. * 0.3 *std::pow(1.-2.872*2.872/(ener*ener),3.)*std::pow(2.872*2.872/(ener*ener),0.8);
961 
962  return sigma;
963  }
964 
966 
967  G4double CrossSectionsStrangeness::NpiToLK(Particle const * const p1, Particle const * const p2) {
968  //
969  // Pion-Nucleon producing Lambda-Kaon cross sections
970  //
971  // ratio
972  // p pi0 -> L K+ (1/2)
973  // p pi- -> L K0 (1)
974 
975 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
976 
977  const Particle *pion;
978  const Particle *nucleon;
980  if(iso == 3 || iso == -3)
981  return 0.;
982 
983  if(p1->isPion()){
984  pion = p1;
985  nucleon = p2;
986  }
987  else{
988  nucleon = p1;
989  pion = p2;
990  }
991  G4double sigma = 0.;
992 
993  if(pion->getType() == PiZero)
994  sigma = 0.5 * p_pimToLK0(pion,nucleon);
995  else
996  sigma = p_pimToLK0(pion,nucleon);
997  return sigma;
998  }
999 
1000  G4double CrossSectionsStrangeness::p_pimToLK0(Particle const * const p1, Particle const * const p2) {
1001 
1002  G4double sigma = 0.;
1003  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1004 
1005  if(pLab < 0.911)
1006  return 0.;
1007 
1008  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);
1009  if(sigma < 0.) return 0;
1010  return sigma;
1011  }
1012 
1013  G4double CrossSectionsStrangeness::NpiToSK(Particle const * const p1, Particle const * const p2) {
1014  //
1015  // Pion-Nucleon producing Sigma-Kaon cross sections
1016  //
1017  // ratio
1018  // p pi+ (5/3) p pi0 (11/6) p pi- (2)
1019  //
1020  // p pi+ -> S+ K+ (10)
1021  // p pi0 -> S+ K0 (6)*
1022  // p pi0 -> S0 K+ (5)
1023  // p pi- -> S0 K0 (6)*
1024  // p pi- -> S- K+ (6)
1025 
1026 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1027 
1028  const Particle *pion;
1029  const Particle *nucleon;
1031 
1032  if(p1->isPion()){
1033  pion = p1;
1034  nucleon = p2;
1035  }
1036  else{
1037  nucleon = p1;
1038  pion = p2;
1039  }
1040  G4double sigma = 0.;
1041 
1042  if(iso == 3 || iso == -3)
1043  sigma = p_pipToSpKp(pion,nucleon);
1044  else if(pion->getType() == PiZero)
1045  sigma = p_pizToSzKp(pion,nucleon)+p_pimToSzKz(pion,nucleon);
1046  else if(iso == 1 || iso == -1)
1047  sigma = p_pimToSzKz(pion,nucleon)+p_pimToSmKp(pion,nucleon);
1048  else // should never append
1049  sigma = 0.;
1050 
1051  return sigma;
1052  }
1053  G4double CrossSectionsStrangeness::p_pimToSmKp(Particle const * const p1, Particle const * const p2) {
1054 
1055  G4double sigma = 0.;
1056  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1057 
1058  if(pLab < 1.0356)
1059  return 0.;
1060 
1061  sigma = 4.352*std::pow(pLab-1.0356,1.006)/(std::pow(pLab+1.0356,0.0978)*std::pow(pLab,5.375));
1062 
1063  if(sigma < 0.) // should never append
1064  return 0;
1065 
1066  return sigma;
1067  }
1068  G4double CrossSectionsStrangeness::p_pipToSpKp(Particle const * const p1, Particle const * const p2) {
1069 
1070  G4double sigma = 0.;
1071  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1072 
1073  if(pLab < 1.0428)
1074  return 0.;
1075 
1076  sigma = 0.001897*std::pow(pLab-1.0428,2.869)/(std::pow(pLab+1.0428,-16.68)*std::pow(pLab,19.1));
1077 
1078  if(sigma < 0.) // should never append
1079  return 0;
1080 
1081  return sigma;
1082  }
1083  G4double CrossSectionsStrangeness::p_pizToSzKp(Particle const * const p1, Particle const * const p2) {
1084 
1085  G4double sigma = 0.;
1086  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1087 
1088  if(pLab < 1.0356)
1089  return 0.;
1090 
1091  sigma = 3.624*std::pow(pLab-1.0356,1.4)/std::pow(pLab,5.14);
1092 
1093  if(sigma < 0.) // should never append
1094  return 0;
1095 
1096  return sigma;
1097  }
1098  G4double CrossSectionsStrangeness::p_pimToSzKz(Particle const * const p1, Particle const * const p2) {
1099 
1100  G4double sigma = 0.;
1101  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(p1,p2); // GeV
1102 
1103  if((p1->getType() == PiZero && pLab < 1.0356) || (pLab < 1.034))
1104  return 0.;
1105 
1106  sigma = 0.3474*std::pow(pLab-1.034,0.07678)/std::pow(pLab,1.627);
1107 
1108  if(sigma < 0.) // should never append
1109  return 0;
1110 
1111  return sigma;
1112  }
1113 
1114  G4double CrossSectionsStrangeness::NpiToLKpi(Particle const * const p1, Particle const * const p2) {
1115  //
1116  // Pion-Nucleon producing Lambda-Kaon-pion cross sections
1117  //
1118  // ratio
1119  // p pi+ (1) p pi0 (3/2) p pi- (2)
1120  //
1121  // p pi0 -> L K+ pi0 (1/2)
1122  // all the others (1)
1123  //
1124 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1125 
1126  G4double sigma=0.;
1127  const Particle *pion;
1128  const Particle *nucleon;
1129 
1130  if(p1->isPion()){
1131  pion = p1;
1132  nucleon = p2;
1133  }
1134  else{
1135  nucleon = p1;
1136  pion = p2;
1137  }
1139  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1140 
1141  if(pLab < 1.147)
1142  return 0.;
1143 
1144  if(iso == 3 || iso == -3)
1145  sigma = 146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1146  else if(pion->getType() == PiZero)
1147  sigma = 1.5*146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1148  else
1149  sigma = 2*146.2*std::pow(pLab-1.147,1.996)/std::pow(pLab+1.147,5.921);
1150 
1151  return sigma;
1152  }
1153 
1154  G4double CrossSectionsStrangeness::NpiToSKpi(Particle const * const p1, Particle const * const p2) {
1155  //
1156  // Pion-Nucleon producing Sigma-Kaon-pion cross sections
1157  //
1158  //ratio
1159  // p pi+ (2.25) p pi0 (2.625) p pi-(3)
1160  //
1161  // p pi+ -> S+ pi+ K0 (5/4)
1162  // p pi+ -> S+ pi0 K+ (3/4)
1163  // p pi+ -> S0 pi+ K+ (1/4)
1164  // p pi0 -> S+ pi0 K0 (1/2)
1165  // p pi0 -> S+ pi- K+ (1/2)
1166  // p pi0 -> S0 pi+ K0 (3/4)
1167  // p pi0 -> S0 pi0 K+ (3/8)
1168  // p pi0 -> S- pi+ K+ (1/2)
1169  // p pi- -> S+ pi- K0 (3/8)
1170  // p pi- -> S0 pi0 K0 (5/8)
1171  // p pi- -> S0 pi- K+ (5/8)
1172  // p pi- -> S- pi+ K0 (1)
1173  // p pi- -> S- pi0 K+ (3/8)
1174 
1175 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1176 
1177  G4double sigma=0.;
1178  const Particle *pion;
1179  const Particle *nucleon;
1180 
1181  if(p1->isPion()){
1182  pion = p1;
1183  nucleon = p2;
1184  }
1185  else{
1186  nucleon = p1;
1187  pion = p2;
1188  }
1190  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1191 
1192  if(pLab <= 1.3041)
1193  return 0.;
1194 
1195  if(iso == 3 || iso == -3)
1196  sigma = 2.25*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1197  else if(pion->getType() == PiZero)
1198  sigma = 2.625*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1199  else
1200  sigma = 3.*8.139*std::pow(pLab-1.3041,2.431)/std::pow(pLab,5.298);
1201 
1202  return sigma;
1203  }
1204 
1205  G4double CrossSectionsStrangeness::NpiToLK2pi(Particle const * const p1, Particle const * const p2) {
1206  //
1207  // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1208  //
1209  // p pi+ (2) p pi0 (1.75) p pi- (2.5)
1210  //
1211  // p pi+ -> L K+ pi+ pi0 (1)
1212  // p pi+ -> L K0 pi+ pi+ (1)
1213  // p pi0 -> L K+ pi0 pi0 (1/4)
1214  // p pi0 -> L K+ pi+ pi- (1)
1215  // p pi0 -> L K0 pi+ pi0 (1/2)
1216  // p pi- -> L K+ pi0 pi- (1)
1217  // p pi- -> L K0 pi+ pi- (1)
1218  // p pi- -> L K0 pi0 pi0 (1/2)
1219 
1220 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1221 
1222  G4double sigma=0.;
1223  const Particle *pion;
1224  const Particle *nucleon;
1225 
1226  if(p1->isPion()){
1227  pion = p1;
1228  nucleon = p2;
1229  }
1230  else{
1231  nucleon = p1;
1232  pion = p2;
1233  }
1235  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1236 
1237  if(pLab <= 1.4162)
1238  return 0.;
1239 
1240  if(iso == 3 || iso == -3)
1241  sigma = 2*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1242  else if(pion->getType() == PiZero)
1243  sigma = 1.75*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1244  else
1245  sigma = 2.5*18.77*std::pow(pLab-1.4162,4.597)/std::pow(pLab,6.877);
1246 
1247  return sigma;
1248  }
1249 
1250  G4double CrossSectionsStrangeness::NpiToSK2pi(Particle const * const p1, Particle const * const p2) {
1251  //
1252  // Pion-Nucleon producing Lambda-Kaon-2pion cross sections
1253  //
1254  // ratio
1255  // p pi+ (3.25) p pi0 (3.5) p pi- (3.75)
1256  //
1257  // p pi+ -> S+ K+ pi+ pi- (1)
1258  // p pi+ -> S+ K+ pi0 pi0 (1/4)
1259  // p pi+ -> S0 K+ pi+ pi0 (1/2)
1260  // p pi+ -> S- K+ pi+ pi+ (1/4)
1261  // p pi+ -> S+ K0 pi+ pi0 (1)
1262  // p pi+ -> S0 K0 pi+ pi+ (1/4)
1263  //
1264  // p pi0 -> S+ K+ pi0 pi- (1/2)
1265  // p pi0 -> S0 K+ pi+ pi- (1/2)
1266  // p pi0 -> S0 K+ pi0 pi0 (1/4)
1267  // p pi0 -> S- K+ pi+ pi0 (1/4)
1268  // p pi0 -> S+ K0 pi+ pi- (1)
1269  // p pi0 -> S+ K0 pi0 pi0 (1/4)
1270  // p pi0 -> S0 K0 pi+ pi0 (1/4)
1271  // p pi0 -> S- K0 pi+ pi+ (1/2)
1272  //
1273  // p pi- -> S+ K+ pi- pi- (1/4)
1274  // p pi- -> S0 K+ pi0 pi- (1/2)
1275  // p pi- -> S- K+ pi+ pi- (1/4)
1276  // p pi- -> S- K+ pi0 pi0 (1/4)
1277  // p pi- -> S+ K0 pi0 pi- (1/2)
1278  // p pi- -> S0 K0 pi+ pi- (1)
1279  // p pi- -> S0 K0 pi0 pi0 (1/2)
1280  // p pi- -> S- K0 pi+ pi0 (1/2)
1281 
1282 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1283 
1284  G4double sigma=0.;
1285  const Particle *pion;
1286  const Particle *nucleon;
1287 
1288  if(p1->isPion()){
1289  pion = p1;
1290  nucleon = p2;
1291  }
1292  else{
1293  nucleon = p1;
1294  pion = p2;
1295  }
1297  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion, nucleon); // GeV
1298 
1299  if(pLab <= 1.5851)
1300  return 0.;
1301 
1302  if(iso == 3 || iso == -3)
1303  sigma = 3.25*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1304  else if(pion->getType() == PiZero)
1305  sigma = 3.5*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1306  else
1307  sigma = 3.75*137.6*std::pow(pLab-1.5851,5.856)/std::pow(pLab,9.295);
1308 
1309  return sigma;
1310  }
1311 
1312  G4double CrossSectionsStrangeness::NpiToNKKb(Particle const * const p1, Particle const * const p2) {
1313  //
1314  // Pion-Nucleon producing Nucleon-Kaon-antiKaon cross sections
1315  //
1316  // ratio
1317  // p pi+ (1/2) p pi0 (3/2) p pi- (5/2)
1318  //
1319  // p pi+ -> p K+ K0b (1/2)
1320  // p pi0 -> p K0 K0b (1/4)
1321  // p pi0 -> p K+ K- (1/4)
1322  // p pi0 -> n K+ K0b (1)
1323  // p pi- -> p K0 K- (1/2)
1324  // p pi- -> n K+ K- (1)
1325  // p pi- -> n K0 K0b (1)
1326 
1327 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1328 
1329  const Particle *particle1;
1330  const Particle *particle2;
1331 
1332  if(p1->isPion()){
1333  particle1 = p1;
1334  particle2 = p2;
1335  }
1336  else{
1337  particle1 = p2;
1338  particle2 = p1;
1339  }
1340 
1341  G4double sigma = 0.;
1342 
1343  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1344 
1345  if(particle1->getType() == PiZero){
1346  if(pLab < 1.5066) return 0.;
1347  else if(pLab < 30.) sigma = 3./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1348  else return 0.;
1349  }
1350  else if((particle1->getType() == PiPlus && particle2->getType() == Neutron) || (particle1->getType() == PiMinus && particle2->getType() == Proton)){
1351  if(pLab < 1.5066) return 0.;
1352  else if(pLab < 30.) sigma = 5./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1353  else return 0.;
1354  }
1355  else{
1356  if(pLab < 1.5066) return 0.;
1357  else if(pLab < 30.) sigma = 1./2.*2.996*std::pow((pLab - 1.5066),1.929)/std::pow(pLab,3.582);
1358  else return 0.;
1359  }
1360  return sigma;
1361  }
1362 
1364  //
1365  // Pion-Nucleon missing strangeness production cross sections
1366  //
1367 // assert((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon()));
1368 
1369  const Particle *pion;
1370  const Particle *nucleon;
1371 
1372  if(p1->isPion()){
1373  pion = p1;
1374  nucleon = p2;
1375  }
1376  else{
1377  pion = p2;
1378  nucleon = p1;
1379  }
1380 
1381  G4double sigma = 0.;
1382 
1383  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(pion,nucleon); // GeV
1384  if(pLab < 2.2) return 0.;
1385 
1386  if(pion->getType() == PiZero){
1387  if(pLab < 30.) sigma = 4.4755*std::pow((pLab - 2.2),1.927)/std::pow(pLab,1.89343);
1388  else return 0.;
1389  }
1390  else if((pion->getType() == PiPlus && nucleon->getType() == Neutron) || (pion->getType() == PiMinus && nucleon->getType() == Proton)){
1391  if(pLab < 30.) sigma = 5.1*std::pow((pLab - 2.2),1.854)/std::pow(pLab,1.904);
1392  else return 0.;
1393  }
1394  else{
1395  if(pLab < 30.) sigma = 3.851*std::pow((pLab - 2.2),2)/std::pow(pLab,1.88286);
1396  else return 0.;
1397  }
1398  return sigma;
1399  }
1400 
1402 
1403  G4double CrossSectionsStrangeness::NLToNS(Particle const * const p1, Particle const * const p2) {
1404  //
1405  // Nucleon-Lambda producing Nucleon-Sigma cross sections
1406  //
1407  // ratio
1408  // p L -> p S0 (1/2)
1409  // p L -> n S+ (1)
1410 
1411 
1412 // assert((p1->isLambda() && p2->isNucleon()) || (p2->isLambda() && p1->isNucleon()));
1413 
1414  G4double sigma = 0.;
1415 
1416  const Particle *particle1;
1417  const Particle *particle2;
1418 
1419  if(p1->isLambda()){
1420  particle1 = p1;
1421  particle2 = p2;
1422  }
1423  else{
1424  particle1 = p2;
1425  particle2 = p1;
1426  }
1427 
1428  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2);
1429 
1430  if(pLab < 0.664)
1431  return 0.;
1432 
1433  sigma = 3 * 8.74*std::pow((pLab-0.664),0.438)/std::pow(pLab,2.717); // 3 * L p -> S0 p
1434 
1435  return sigma;
1436  }
1437 
1438  G4double CrossSectionsStrangeness::NSToNL(Particle const * const p1, Particle const * const p2) {
1439  //
1440  // Nucleon-Lambda producing Nucleon-Sigma cross sections
1441  //
1442  // ratio
1443  // p S0 -> p L (1/2)
1444  // p S- -> n L (1)
1445 
1446 // assert((p1->isSigma() && p2->isNucleon()) || (p2->isSigma() && p1->isNucleon()));
1447 
1449  if(iso == 3 || iso == -3)
1450  return 0.;
1451 
1452  G4double sigma;
1453  const Particle *particle1;
1454  const Particle *particle2;
1455 
1456  if(p1->isSigma()){
1457  particle1 = p1;
1458  particle2 = p2;
1459  }
1460  else{
1461  particle1 = p2;
1462  particle2 = p1;
1463  }
1464  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1465 
1466  if(particle1->getType() == SigmaZero){
1467  if(pLab < 0.1) return 100.; // cut-max
1468  sigma = 8.23*std::pow(pLab,-1.087);
1469  }
1470  else{
1471  if(pLab < 0.1) return 200.; // cut-max
1472  sigma = 16.46*std::pow(pLab,-1.087);
1473  }
1474  return sigma;
1475  }
1476 
1477  G4double CrossSectionsStrangeness::NSToNS(Particle const * const p1, Particle const * const p2) {
1478 
1479 // assert((p1->isSigma() && p2->isNucleon()) || (p2->isSigma() && p1->isNucleon()));
1480 
1482  if(iso == 3 || iso == -3)
1483  return 0.; // only quasi-elastic here
1484 
1485  const Particle *particle1;
1486  const Particle *particle2;
1487 
1488  if(p1->isSigma()){
1489  particle1 = p1;
1490  particle2 = p2;
1491  }
1492  else{
1493  particle1 = p2;
1494  particle2 = p1;
1495  }
1496  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1497 
1498  if(particle2->getType() == Neutron && pLab < 0.162) return 0.;
1499  else if(pLab < 0.1035) return 200.; // cut-max
1500 
1501  return 13.79*std::pow(pLab,-1.181);
1502  }
1503 
1505 
1506  G4double CrossSectionsStrangeness::NKToNK(Particle const * const p1, Particle const * const p2) {
1507  //
1508  // Nucleon-Kaon quasi-elastic cross sections
1509  //
1510 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1511 
1513 
1514  if(iso != 0)
1515  return 0.;
1516 
1517  const Particle *particle1;
1518  const Particle *particle2;
1519 
1520  if(p1->isKaon()){
1521  particle1 = p1;
1522  particle2 = p2;
1523  }
1524  else{
1525  particle1 = p2;
1526  particle2 = p1;
1527  }
1528 
1529  G4double sigma = 0.;
1530  G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1531 
1532  if(particle1->getType() == Proton)
1533  pLab += 2*0.0774;
1534 
1535  if(pLab <= 0.0774)
1536  return 0.;
1537 
1538  sigma = 12.84*std::pow((pLab-0.0774),18.19)/std::pow((pLab),20.41);
1539 
1540  return sigma;
1541  }
1542 
1543  G4double CrossSectionsStrangeness::NKToNKpi(Particle const * const p1, Particle const * const p2) {
1544  //
1545  // Nucleon-Kaon producing Nucleon-Kaon-pion cross sections
1546  //
1547  // Ratio determined by meson symmetry using only "resonante" diagram (with Delta or K*)
1548  //
1549  // ratio: K+ p (5) K0 p (5.545)
1550  //
1551  // K+ p -> p K+ pi0 1.2
1552  // K+ p -> p K0 pi+ 3
1553  // K+ p -> n K+ pi+ 0.8
1554  // K0 p -> p K+ pi- 1
1555  // K0 p -> p K0 pi0 0.845
1556  // K0 p -> n K+ pi0 1.47
1557  // K0 p -> n K0 pi+ 2.23
1558 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1559 
1561 
1562  const Particle *particle1;
1563  const Particle *particle2;
1564 
1565  if(p1->isKaon()){
1566  particle1 = p1;
1567  particle2 = p2;
1568  }
1569  else{
1570  particle1 = p2;
1571  particle2 = p1;
1572  }
1573 
1574  G4double sigma = 0.;
1575  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1576 
1577  if(pLab <= 0.53)
1578  return 0.;
1579 
1580  if(iso == 0)
1581  sigma = 5.55*116.8*std::pow(pLab-0.53,6.874)/std::pow(pLab,10.11);
1582  else
1583  sigma = 5.*116.8*std::pow(pLab-0.53,6.874)/std::pow(pLab,10.11);;
1584  return sigma;
1585  }
1586 
1587  G4double CrossSectionsStrangeness::NKToNK2pi(Particle const * const p1, Particle const * const p2) {
1588  //
1589  // Nucleon-Kaon producing Nucleon-Kaon-2pion cross sections
1590  //
1591  // p K+ (2.875) p K0 (3.125)
1592  //
1593  // p K+ -> p K+ pi+ pi- (1)
1594  // p K+ -> p K+ pi0 pi0 (1/8)
1595  // p K+ -> p K0 pi+ pi0 (1)
1596  // p K+ -> n K+ pi+ pi0 (1/2)
1597  // p K+ -> n K0 pi+ pi+ (1/4)
1598  // p K0 -> p K+ pi0 pi- (1)
1599  // p K0 -> p K0 pi+ pi- (1)
1600  // p K0 -> p K0 pi0 pi0 (1/8)
1601  // p K0 -> n K+ pi+ pi- (1/4)
1602  // p K0 -> n K+ pi0 pi0 (1/4)
1603  // p K0 -> n K0 pi+ pi0 (1/2)
1604 
1605 // assert((p1->isNucleon() && p2->isKaon()) || (p2->isNucleon() && p1->isKaon()));
1606 
1608 
1609  const Particle *particle1;
1610  const Particle *particle2;
1611 
1612  if(p1->isKaon()){
1613  particle1 = p1;
1614  particle2 = p2;
1615  }
1616  else{
1617  particle1 = p2;
1618  particle2 = p1;
1619  }
1620 
1621  G4double sigma = 0.;
1622  const G4double pLab = 0.001 * KinematicsUtils::momentumInLab(particle1,particle2); // GeV
1623 
1624  if(pLab < 0.812)
1625  sigma = 0.;
1626  else if(pLab < 1.744)
1627  sigma = 26.41*std::pow(pLab-0.812,7.138)/std::pow(pLab,5.337);
1628  else if(pLab < 3.728)
1629  sigma = 1572.*std::pow(pLab-0.812,9.069)/std::pow(pLab,12.44);
1630  else
1631  sigma = 60.23*std::pow(pLab-0.812,5.084)/std::pow(pLab,6.72);
1632 
1633  if(iso == 0)
1634  sigma *= 3.125;
1635  else
1636  sigma *= 2.875;
1637 
1638  return sigma;
1639  }
1640 
1642 
1643  G4double CrossSectionsStrangeness::NKbToNKb(Particle const * const p1, Particle const * const p2) {
1644  //
1645  // Nucleon-antiKaon quasi-elastic cross sections
1646  //
1647 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1648 
1649  G4double sigma=0.;
1651 
1652  const Particle *antikaon;
1653  const Particle *nucleon;
1654 
1655  if (p1->isAntiKaon()) {
1656  antikaon = p1;
1657  nucleon = p2;
1658  }
1659  else {
1660  antikaon = p2;
1661  nucleon = p1;
1662  }
1663 
1664  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1665 
1666  if(iso != 0) // K0b p and K- n -> forbidden: quasi-elastic diffusion only
1667  return 0;
1668  else if(nucleon->getType() == Proton){ // K- p -> K0b n
1669  if(pLab < 0.08921)
1670  return 0.;
1671  else if(pLab < 0.2)
1672  sigma = 0.4977*std::pow(pLab - 0.08921,0.5581)/std::pow(pLab,2.704);
1673  else if(pLab < 0.73)
1674  sigma = 2.*std::pow(pLab,-1.2) + 6.493*std::exp(-0.5*std::pow((pLab-0.3962)/0.02,2));
1675  else if(pLab < 1.38)
1676  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));
1677  else if(pLab < 30.)
1678  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));
1679  else sigma = 0.;
1680  }
1681  else{ // K0b n -> K- p (same as K- p but without threshold)
1682  if(pLab < 0.1)
1683  sigma = 30.;
1684  else if(pLab < 0.73)
1685  sigma = 2.*std::pow(pLab,-1.2) + 6.493*std::exp(-0.5*std::pow((pLab-0.3962)/0.02,2));
1686  else if(pLab < 1.38)
1687  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));
1688  else if(pLab < 30.)
1689  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));
1690  else sigma = 0.;
1691  }
1692  return sigma;
1693  }
1694 
1695  G4double CrossSectionsStrangeness::NKbToSpi(Particle const * const p1, Particle const * const p2) {
1696  //
1697  // Nucleon-antiKaon producing Sigma-pion cross sections
1698  //
1699  // ratio
1700  // p K0b (4/3) p K- (13/6)
1701  //
1702  // p K0b -> S+ pi0 (2/3)
1703  // p K0b -> S0 pi+ (2/3)
1704  // p K- -> S+ pi- (1)
1705  // p K- -> S0 pi0 (1/2)
1706  // p K- -> S- pi+ (2/3)
1707 
1708 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1709 
1710  G4double sigma=0.;
1712 
1713  const Particle *antikaon;
1714  const Particle *nucleon;
1715 
1716  if (p1->isAntiKaon()) {
1717  antikaon = p1;
1718  nucleon = p2;
1719  }
1720  else {
1721  antikaon = p2;
1722  nucleon = p1;
1723  }
1724 
1725  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1726 
1727  if(iso == 0){
1728  if(pLab < 0.1)
1729  return 152.0; // 70.166*13/6
1730  else
1731  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));
1732  }
1733  else{
1734  if(pLab < 0.1)
1735  return 93.555; // 70.166*4/3
1736  else
1737  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));
1738  //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));
1739  }
1740 
1741  return sigma;
1742  }
1743 
1744  G4double CrossSectionsStrangeness::NKbToLpi(Particle const * const p1, Particle const * const p2) {
1745  //
1746  // Nucleon-antiKaon producing Lambda-pion cross sections
1747  //
1748  // ratio
1749  // p K0b (1) p K- (1/2)
1750  //
1751  // p K- -> L pi0 (1/2)
1752  // p K0b -> L pi+ (1)
1753 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1754 
1755  G4double sigma = 0.;
1757 
1758  const Particle *antikaon;
1759  const Particle *nucleon;
1760 
1761  if (p1->isAntiKaon()) {
1762  antikaon = p1;
1763  nucleon = p2;
1764  }
1765  else {
1766  antikaon = p2;
1767  nucleon = p1;
1768  }
1769  if(iso == 0)
1770  sigma = p_kmToL_pz(antikaon,nucleon);
1771  else
1772  sigma = 2*p_kmToL_pz(antikaon,nucleon);
1773 
1774  return sigma;
1775  }
1776  G4double CrossSectionsStrangeness::p_kmToL_pz(Particle const * const p1, Particle const * const p2) {
1777  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(p1, p2); // GeV
1778  G4double sigma = 0.;
1779  if(pLab < 0.086636)
1780  sigma = 40.24;
1781  else if(pLab < 0.5)
1782  sigma = 0.97*std::pow(pLab,-1.523);
1783  else if(pLab < 2.)
1784  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);
1785  else if(pLab < 30.)
1786  sigma = 3.*std::pow(pLab,-2.57);
1787  else
1788  sigma = 0.;
1789 
1790  return sigma;
1791  }
1792 
1793  G4double CrossSectionsStrangeness::NKbToS2pi(Particle const * const p1, Particle const * const p2) {
1794  //
1795  // Nucleon-antiKaon producing Sigma-2pion cross sections
1796  //
1797  // ratio
1798  // p K0b (29/12) p K- (59/24)
1799  //
1800  // p K0b -> S+ pi+ pi- (2/3)
1801  // p K0b -> S+ pi0 pi0 (1/4)
1802  // p K0b -> S0 pi+ pi0 (5/6)
1803  // p K0b -> S- pi+ pi+ (2/3)
1804  // p K- -> S+ pi0 pi- (1)
1805  // p K- -> S0 pi+ pi- (2/3)
1806  // p K- -> S0 pi0 pi0 (1/8)
1807  // p K- -> S- pi+ pi0 (2/3)
1808 
1809 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1810 
1811  G4double sigma=0.;
1813 
1814  const Particle *antikaon;
1815  const Particle *nucleon;
1816 
1817  if (p1->isAntiKaon()) {
1818  antikaon = p1;
1819  nucleon = p2;
1820  }
1821  else {
1822  antikaon = p2;
1823  nucleon = p1;
1824  }
1825 
1826  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1827 
1828  if(pLab < 0.260)
1829  return 0.;
1830 
1831  if(iso == 0)
1832  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));
1833  else
1834  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));
1835 
1836  /*
1837  if(iso == 0)
1838  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));
1839  else
1840  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));*/
1841 
1842  return sigma;
1843  }
1844 
1845  G4double CrossSectionsStrangeness::NKbToL2pi(Particle const * const p1, Particle const * const p2) {
1846  //
1847  // Nucleon-antiKaon producing Lambda-2pion cross sections
1848  //
1849  // ratio
1850  // p K0b -> L pi+ pi0 (1)
1851  // p K- -> L pi+ pi- (1)
1852  // p K- -> L pi0 pi0 (1/4)
1853 
1854 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1855 
1856  G4double sigma = 0.;
1858 
1859  const Particle *antikaon;
1860  const Particle *nucleon;
1861 
1862  if (p1->isAntiKaon()) {
1863  antikaon = p1;
1864  nucleon = p2;
1865  }
1866  else {
1867  antikaon = p2;
1868  nucleon = p1;
1869  }
1870 
1871  if(iso == 0)
1872  sigma = 1.25*p_kmToL_pp_pm(antikaon,nucleon);
1873  else
1874  sigma = p_kmToL_pp_pm(antikaon,nucleon);
1875 
1876  return sigma;
1877  }
1879  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(p1, p2); // GeV
1880  G4double sigma = 0.;
1881  if(pLab < 0.97)
1882  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.);
1883  else if(pLab < 30)
1884  sigma = 46.3*std::pow(pLab,0.62)/std::pow(pLab+1.,3.565);
1885  else
1886  sigma = 0.;
1887 
1888  return sigma;
1889  }
1890 
1891  G4double CrossSectionsStrangeness::NKbToNKbpi(Particle const * const p1, Particle const * const p2) {
1892  //
1893  // Nucleon-antiKaon producing Nucleon-antiKaon-pion cross sections
1894  //
1895  // ratio
1896  // p K0b (2) p K- (4)
1897  //
1898  // p K0b -> p K0b pi0 (1/2)
1899  // p K0b -> p K- pi+ (1)*
1900  // p K0b -> n K0b pi+ (1/2)*
1901  // p K- -> p K- pi0 (1/2)*
1902  // p K- -> p K0b pi- (2/3)*
1903  // p K- -> n K- pi+ (3/4)*
1904  // p K- -> n K0b pi0 (2)
1905 
1906  //
1907 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1908 
1909  G4double sigma=0.;
1911 
1912  const Particle *antikaon;
1913  const Particle *nucleon;
1914 
1915  if (p1->isAntiKaon()) {
1916  antikaon = p1;
1917  nucleon = p2;
1918  }
1919  else {
1920  antikaon = p2;
1921  nucleon = p1;
1922  }
1923 
1924  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1925 
1926  if(pLab < 0.526)
1927  return 0.;
1928 
1929  if(iso == 0)
1930  sigma = 2. * 101.3*std::pow(pLab-0.526,5.846)/std::pow(pLab,8.343);
1931  else
1932  sigma = 4. * 101.3*std::pow(pLab-0.526,5.846)/std::pow(pLab,8.343);
1933 
1934  return sigma;
1935  }
1936 
1937  G4double CrossSectionsStrangeness::NKbToNKb2pi(Particle const * const p1, Particle const * const p2) {
1938  //
1939  // Nucleon-antiKaon producing Nucleon-antiKaon-2pion cross sections
1940  //
1941  // ratio
1942  // p K0b (4.25) p K- (4.75)
1943  //
1944  // p K0b -> p K0b pi+ pi- (1)
1945  // p K0b -> p K0b pi0 pi0 (1/4)
1946  // p K0b -> p K- pi+ pi0 (1)
1947  // p K0b -> n K0b pi+ pi0 (1)
1948  // p K0b -> n K- pi+ pi+ (1)
1949  // p K- -> p K0b pi0 pi- (1)
1950  // p K- -> p K- pi+ pi- (1)
1951  // p K- -> p K- pi0 pi0 (1/4)
1952  // p K- -> n K0b pi+ pi- (1)
1953  // p K- -> n K0b pi0 pi0 (1/2)
1954  // p K- -> n K- pi+ pi0 (1)
1955 
1956 // assert((p1->isNucleon() && p2->isAntiKaon()) || (p1->isAntiKaon() && p2->isNucleon()));
1957 
1958  G4double sigma=0.;
1960 
1961  const Particle *antikaon;
1962  const Particle *nucleon;
1963 
1964  if (p1->isAntiKaon()) {
1965  antikaon = p1;
1966  nucleon = p2;
1967  }
1968  else {
1969  antikaon = p2;
1970  nucleon = p1;
1971  }
1972 
1973  const G4double pLab = 0.001*KinematicsUtils::momentumInLab(antikaon, nucleon); // GeV
1974 
1975  if(pLab < 0.85)
1976  return 0.;
1977 
1978  if(iso == 0)
1979  sigma = 4.75 * 26.8*std::pow(pLab-0.85,4.9)/std::pow(pLab,6.34);
1980  else
1981  sigma = 4.25 * 26.8*std::pow(pLab-0.85,4.9)/std::pow(pLab,6.34);
1982 
1983  return sigma;
1984  }
1985 
1986 
1987 } // namespace G4INCL
1988 
G4double NNTot(Particle const *const part1, Particle const *const part2)
Internal implementation of the NN total cross section.
Float_t x
Definition: compare.C:6
virtual G4double NNToNSK(Particle const *const p1, Particle const *const p2)
G4bool pion(G4int ityp)
virtual G4double NKbToL2pi(Particle const *const p1, Particle const *const p2)
virtual G4double elastic(Particle const *const p1, Particle const *const p2)
second new elastic particle-particle cross section
virtual G4double NNToxPiNN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NN Channel.
virtual G4double p_pimToSzKz(Particle const *const p1, Particle const *const p2)
const G4double effectiveNucleonMass
virtual G4double NNToNSKpi(Particle const *const p1, Particle const *const p2)
virtual G4double p_pimToSmKp(Particle const *const p1, Particle const *const p2)
virtual G4double NpiToLKpi(Particle const *const p1, Particle const *const p2)
virtual G4double NpiToSK2pi(Particle const *const p1, Particle const *const p2)
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
G4bool isOmega() const
Is this an omega?
virtual G4double NpiToSKpi(Particle const *const p1, Particle const *const p2)
virtual G4double NKbToS2pi(Particle const *const p1, Particle const *const p2)
virtual G4double NNToNLK2pi(Particle const *const p1, Particle const *const p2)
virtual G4double NDeltaToNNKKb(Particle const *const p1, Particle const *const p2)
virtual G4double NSToNS(Particle const *const p1, Particle const *const p2)
G4double squareTotalEnergyInCM(Particle const *const p1, Particle const *const p2)
G4bool nucleon(G4int ityp)
G4bool isHyperon() const
Is this an Hyperon?
virtual G4double etaNElastic(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - elastic Channel.
virtual G4double NLToNS(Particle const *const p1, Particle const *const p2)
Nucleon-Hyperon quasi-elastic cross sections.
G4double p_pipToSpKp(Particle const *const p1, Particle const *const p2)
G4double getMass() const
Get the cached particle mass.
G4bool isSigma() const
Is this a Sigma?
virtual G4double NSToNL(Particle const *const p1, Particle const *const p2)
virtual G4double NDeltaToNLK(Particle const *const p1, Particle const *const p2)
Nucleon-Delta to Stange particles cross sections.
virtual G4double NKbelastic(Particle const *const p1, Particle const *const p2)
virtual G4double p_pizToSzKp(Particle const *const p1, Particle const *const p2)
G4double p_pimToLK0(Particle const *const p1, Particle const *const p2)
virtual G4double etaNToPiPiN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - pipiN Channel.
virtual G4double NKToNK2pi(Particle const *const p1, Particle const *const p2)
G4bool isLambda() const
Is this a Lambda?
virtual G4double NNToNNEta(Particle const *const particle1, Particle const *const particle2)
Cross section for Eta production (inclusive) - NN entrance channel.
virtual G4double NKbToNKb2pi(Particle const *const p1, Particle const *const p2)
static const G4int nMaxPiPiN
Maximum number of outgoing pions in piN collisions.
virtual G4double NNToNNKKb(Particle const *const p1, Particle const *const p2)
const G4double effectiveNucleonMass2
const XML_Char * s
Definition: expat.h:262
double G4double
Definition: G4Types.hh:76
virtual G4double NpiToLK(Particle const *const p1, Particle const *const p2)
Nucleon-Pion to Stange particles cross sections.
virtual G4double NYelastic(Particle const *const p1, Particle const *const p2)
elastic scattering for Nucleon-Strange Particles cross sections
virtual G4double NNToNLK(Particle const *const p1, Particle const *const p2)
Nucleon-Nucleon to Stange particles cross sections.
virtual G4double total(Particle const *const p1, Particle const *const p2)
second new total particle-particle cross section
G4bool isDelta() const
Is it a Delta?
virtual G4double etaNToPiN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - piN Channel.
virtual G4double NDeltaToDeltaSK(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 NNToNNOmega(Particle const *const particle1, Particle const *const particle2)
Cross section for Omega production (inclusive) - NN entrance channel.
static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients< N > const &coeffs)
virtual G4double NpiToNKKb(Particle const *const p1, Particle const *const p2)
virtual G4double NKbToLpi(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
virtual G4double elastic(Particle const *const p1, Particle const *const p2)
Elastic particle-particle cross section.
virtual G4double piNToEtaN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance production - piN Channel.
Multipion, mesonic Resonances and strange cross sections.
virtual G4double NDeltaToNSK(Particle const *const p1, Particle const *const p2)
virtual G4double NKToNKpi(Particle const *const p1, Particle const *const p2)
virtual G4double NKbToNKbpi(Particle const *const p1, Particle const *const p2)
virtual G4double NpiToSK(Particle const *const p1, Particle const *const p2)
virtual G4double NKToNK(Particle const *const p1, Particle const *const p2)
Nucleon-Kaon cross sections.
G4double piNTot(Particle const *const p1, Particle const *const p2)
G4bool hyperon(G4int ityp)
G4double NNInelasticIso(const G4double ener, const G4int iso)
Internal implementation of the isospin dependent NN reaction cross section.
virtual G4double NpiToLK2pi(Particle const *const p1, Particle const *const p2)
virtual G4double NKbToNKb(Particle const *const p1, Particle const *const p2)
Nucleon-antiKaon cross sections.
static G4double eval(const G4double x, HornerCoefficients< N > const &coeffs)
int G4int
Definition: G4Types.hh:78
virtual G4double etaPrimeNToPiN(Particle const *const p1, Particle const *const p2)
Cross section for EtaPrimeN-&gt;PiN.
virtual G4double NKbToSpi(Particle const *const p1, Particle const *const p2)
static const G4int nMaxPiNN
Maximum number of outgoing pions in NN collisions.
virtual G4double NNTwoPi(const G4double ener, const G4int iso, const G4double xsiso)
Cross section for direct 2-pion production - NN entrance channel.
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
G4INCL::ParticleType getType() const
G4bool isKaon() const
Is this a Kaon?
virtual G4double NNToNSK2pi(Particle const *const p1, Particle const *const p2)
G4bool isEta() const
Is this an eta?
virtual G4double p_kmToL_pz(Particle const *const p1, Particle const *const p2)
G4bool isPion() const
Is this a pion?
G4bool isAntiKaon() const
Is this an antiKaon?
virtual G4double NpiToMissingStrangeness(Particle const *const p1, Particle const *const p2)
virtual G4double omegaNInelastic(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - inelastic Channel.
virtual G4double NDeltaToDeltaLK(Particle const *const p1, Particle const *const p2)
virtual G4double p_kmToL_pp_pm(Particle const *const p1, Particle const *const p2)
virtual G4double piNToOmegaN(Particle const *const p1, Particle const *const p2)
Cross section for PiN-&gt;OmegaN.
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 NKelastic(Particle const *const p1, Particle const *const p2)
virtual G4double NNToNLKpi(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.
virtual G4double NNToMissingStrangeness(Particle const *const p1, Particle const *const p2)
G4bool isNucleon() const
virtual G4double piNToxPiN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - piN Channel.
virtual G4double NDeltaToNN(Particle const *const p1, Particle const *const p2)
Cross section for NDelta-&gt;NN.
G4bool isEtaPrime() const
Is this an etaprime?