source: trunk/source/processes/hadronic/cross_sections/src/G4HadronNucleonXsc.cc@ 1201

Last change on this file since 1201 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 38.6 KB
RevLine 
[819]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//
27// 14.03.07 V. Grichine - first implementation
28//
29
30#include "G4HadronNucleonXsc.hh"
31
32#include "G4ParticleTable.hh"
33#include "G4IonTable.hh"
34#include "G4ParticleDefinition.hh"
35
36//////////////////////////////////////////////////////////////////////////////////////
37//
38//
39
40
41G4HadronNucleonXsc::G4HadronNucleonXsc()
42: fUpperLimit( 10000 * GeV ),
43 fLowerLimit( 0.03 * MeV )
44{
45 theGamma = G4Gamma::Gamma();
46 theProton = G4Proton::Proton();
47 theNeutron = G4Neutron::Neutron();
48 theAProton = G4AntiProton::AntiProton();
49 theANeutron = G4AntiNeutron::AntiNeutron();
50 thePiPlus = G4PionPlus::PionPlus();
51 thePiMinus = G4PionMinus::PionMinus();
52 thePiZero = G4PionZero::PionZero();
53 theKPlus = G4KaonPlus::KaonPlus();
54 theKMinus = G4KaonMinus::KaonMinus();
55 theK0S = G4KaonZeroShort::KaonZeroShort();
56 theK0L = G4KaonZeroLong::KaonZeroLong();
57 theL = G4Lambda::Lambda();
58 theAntiL = G4AntiLambda::AntiLambda();
59 theSPlus = G4SigmaPlus::SigmaPlus();
60 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
61 theSMinus = G4SigmaMinus::SigmaMinus();
62 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
63 theS0 = G4SigmaZero::SigmaZero();
64 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
65 theXiMinus = G4XiMinus::XiMinus();
66 theXi0 = G4XiZero::XiZero();
67 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
68 theAXi0 = G4AntiXiZero::AntiXiZero();
69 theOmega = G4OmegaMinus::OmegaMinus();
70 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
71 theD = G4Deuteron::Deuteron();
72 theT = G4Triton::Triton();
73 theA = G4Alpha::Alpha();
74 theHe3 = G4He3::He3();
75}
76
77///////////////////////////////////////////////////////////////////////////////////////
78//
79//
80
81G4HadronNucleonXsc::~G4HadronNucleonXsc()
82{
83}
84
85
86////////////////////////////////////////////////////////////////////////////////////////
87//
88//
89
90
91G4bool
92G4HadronNucleonXsc::IsApplicable(const G4DynamicParticle* aDP,
93 const G4Element* anElement)
94{
95 return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN());
96}
97
98////////////////////////////////////////////////////////////////////////////////////////
99//
100//
101
102G4bool
103G4HadronNucleonXsc::IsZAApplicable(const G4DynamicParticle* aDP,
104 G4double Z, G4double)
105{
106 G4bool applicable = false;
107 // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
108 G4double kineticEnergy = aDP->GetKineticEnergy();
109
110 const G4ParticleDefinition* theParticle = aDP->GetDefinition();
111
112 if ( ( kineticEnergy >= fLowerLimit &&
113 Z > 1.5 && // >= He
114 ( theParticle == theAProton ||
115 theParticle == theGamma ||
116 theParticle == theKPlus ||
117 theParticle == theKMinus ||
118 theParticle == theSMinus) ) ||
119
120 ( kineticEnergy >= 0.1*fLowerLimit &&
121 Z > 1.5 && // >= He
122 ( theParticle == theProton ||
123 theParticle == theNeutron ||
124 theParticle == thePiPlus ||
125 theParticle == thePiMinus ) ) ) applicable = true;
126
127 return applicable;
128}
129
130
131
132
133/////////////////////////////////////////////////////////////////////////////////////
134//
135// Returns hadron-nucleon Xsc according to differnt parametrisations:
136// [2] E. Levin, hep-ph/9710546
137// [3] U. Dersch, et al, hep-ex/9910052
138// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
139
140G4double
141G4HadronNucleonXsc::GetHadronNucleonXscEL(const G4DynamicParticle* aParticle,
142 const G4ParticleDefinition* nucleon )
143{
144 G4double xsection;
145
146
147 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
148
149 G4double proj_mass = aParticle->GetMass();
150 G4double proj_momentum = aParticle->GetMomentum().mag();
151 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
152
153 sMand /= GeV*GeV; // in GeV for parametrisation
154 proj_momentum /= GeV;
155
156 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
157
158 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
159
160
161 if(theParticle == theGamma && pORn )
162 {
163 xsection = (0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
164 }
165 else if(theParticle == theNeutron && pORn ) // as proton ???
166 {
167 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
168 }
169 else if(theParticle == theProton && pORn )
170 {
171 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
172
173 // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
174 // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
175 }
176 else if(theParticle == theAProton && pORn )
177 {
178 xsection = ( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
179 }
180 else if(theParticle == thePiPlus && pORn )
181 {
182 xsection = (13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
183 }
184 else if(theParticle == thePiMinus && pORn )
185 {
186 // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
187 xsection = (13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
188 }
189 else if(theParticle == theKPlus && pORn )
190 {
191 xsection = (11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
192 }
193 else if(theParticle == theKMinus && pORn )
194 {
195 xsection = (11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
196 }
197 else // as proton ???
198 {
199 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
200 }
201 xsection *= millibarn;
202
203 fTotalXsc = xsection;
204 fInelasticXsc = 0.83*xsection;
205 fElasticXsc = fTotalXsc - fInelasticXsc;
206 if (fElasticXsc < 0.)fElasticXsc = 0.;
207
208 return xsection;
209}
210
211
212
213
214
215/////////////////////////////////////////////////////////////////////////////////////
216//
217// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
218// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
219// At = number of nucleons, Zt = number of protons
220
221G4double
222G4HadronNucleonXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
223 const G4ParticleDefinition* nucleon )
224{
225 G4double xsection(0);
226 G4int Zt=1, Nt=1, At=1;
227
228 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
229
230 G4double proj_mass = aParticle->GetMass();
231 G4double proj_momentum = aParticle->GetMomentum().mag();
232
233 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
234
235 sMand /= GeV*GeV; // in GeV for parametrisation
236
237 // General PDG fit constants
238
239 G4double s0 = 5.38*5.38; // in Gev^2
240 G4double eta1 = 0.458;
241 G4double eta2 = 0.458;
242 G4double B = 0.308;
243
244
245 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
246
247 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
248 G4bool proton = (nucleon == theProton);
249 G4bool neutron = (nucleon == theNeutron);
250
251
252 if(theParticle == theNeutron) // proton-neutron fit
253 {
254 if ( proton )
255 {
256 xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
257 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));// on p
258 }
259 if ( neutron )
260 {
261 xsection = Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
262 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // on n pp for nn
263 }
264 }
265 else if(theParticle == theProton)
266 {
267 if ( proton )
268 {
269 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
270 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
271 }
272 if ( neutron )
273 {
274 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
275 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
276 }
277 }
278 else if(theParticle == theAProton)
279 {
280 if ( proton )
281 {
282 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
283 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
284 }
285 if ( neutron )
286 {
287 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
288 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
289 }
290 }
291 else if(theParticle == thePiPlus && pORn )
292 {
293 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
294 + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
295 }
296 else if(theParticle == thePiMinus && pORn )
297 {
298 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
299 + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
300 }
301 else if(theParticle == theKPlus)
302 {
303 if ( proton )
304 {
305 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
306 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
307 }
308 if ( neutron )
309 {
310 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
311 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
312 }
313 }
314 else if(theParticle == theKMinus)
315 {
316 if ( proton )
317 {
318 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
319 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
320 }
321 if ( neutron )
322 {
323 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
324 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2) );
325 }
326 }
327 else if(theParticle == theSMinus && pORn )
328 {
329 xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
330 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2) );
331 }
332 else if(theParticle == theGamma && pORn ) // modify later on
333 {
334 xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
335 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2) );
336
337 }
338 else // as proton ???
339 {
340 if ( proton )
341 {
342 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
343 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2) );
344 }
345 if ( neutron )
346 {
347 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
348 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
349 }
350 }
351 xsection *= millibarn; // parametrised in mb
352
353 fTotalXsc = xsection;
354 fInelasticXsc = 0.83*xsection;
355 fElasticXsc = fTotalXsc - fInelasticXsc;
356 if (fElasticXsc < 0.)fElasticXsc = 0.;
357
358 return xsection;
359}
360
361
362
363/////////////////////////////////////////////////////////////////////////////////////
364//
365// Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
366// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
367
368G4double
369G4HadronNucleonXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
370 const G4ParticleDefinition* nucleon )
371{
372 G4double xsection(0), Delta, A0, B0;
373 G4int Zt=1, Nt=1, At=1;
374 G4double hpXsc(0);
375 G4double hnXsc(0);
376
377
378 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
379
380 G4double proj_mass = aParticle->GetMass();
381 G4double proj_energy = aParticle->GetTotalEnergy();
382 G4double proj_momentum = aParticle->GetMomentum().mag();
383
384 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
385
386 sMand /= GeV*GeV; // in GeV for parametrisation
387 proj_momentum /= GeV;
388 proj_energy /= GeV;
389 proj_mass /= GeV;
390
391 // General PDG fit constants
392
393 G4double s0 = 5.38*5.38; // in Gev^2
394 G4double eta1 = 0.458;
395 G4double eta2 = 0.458;
396 G4double B = 0.308;
397
398
399 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
400
401 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
402 G4bool proton = (nucleon == theProton);
403 G4bool neutron = (nucleon == theNeutron);
404
405 if( theParticle == theNeutron && pORn)
406 {
407 if( proj_momentum >= 10.)
408 // if( proj_momentum >= 2.)
409 {
410 Delta = 1.;
411
412 if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
413
414 if(proj_momentum >= 10.)
415 {
416 B0 = 7.5;
417 A0 = 100. - B0*std::log(3.0e7);
418
419 xsection = A0 + B0*std::log(proj_energy) - 11
420 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
421 0.93827*0.93827,-0.165); // mb
422 }
423 fTotalXsc = xsection;
424 }
425 else
426 {
427 // nn to be pp
428
429 if(neutron)
430 {
431 if( proj_momentum < 0.73 )
432 {
433 hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
434 }
435 else if( proj_momentum < 1.05 )
436 {
437 hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
438 (std::log(proj_momentum/0.73));
439 }
440 else // if( proj_momentum < 10. )
441 {
442 hnXsc = 39.0+
443 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
444 }
445 fTotalXsc = hnXsc;
446 }
447 // pn to be np
448
449 if(proton)
450 {
451 if( proj_momentum < 0.8 )
452 {
453 hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
454 }
455 else if( proj_momentum < 1.4 )
456 {
457 hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
458 }
459 else // if( proj_momentum < 10. )
460 {
461 hpXsc = 33.3+
462 20.8*(std::pow(proj_momentum,2.0)-1.35)/
463 (std::pow(proj_momentum,2.50)+0.95);
464 }
465 fTotalXsc = hpXsc;
466 }
467 // xsection = hpXsc*Zt + hnXsc*Nt;
468 }
469 }
470 else if(theParticle == theProton && pORn)
471 {
472 if( proj_momentum >= 10.)
473 // if( proj_momentum >= 2.)
474 {
475 Delta = 1.;
476
477 if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
478
479 if(proj_momentum >= 10.)
480 {
481 B0 = 7.5;
482 A0 = 100. - B0*std::log(3.0e7);
483
484 xsection = A0 + B0*std::log(proj_energy) - 11
485 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
486 0.93827*0.93827,-0.165); // mb
487 }
488 fTotalXsc = xsection;
489 }
490 else
491 {
492 // pp
493
494 if(proton)
495 {
496 if( proj_momentum < 0.73 )
497 {
498 hpXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
499 }
500 else if( proj_momentum < 1.05 )
501 {
502 hpXsc = 23 + 40*(std::log(proj_momentum/0.73))*
503 (std::log(proj_momentum/0.73));
504 }
505 else // if( proj_momentum < 10. )
506 {
507 hpXsc = 39.0+
508 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
509 }
510 fTotalXsc = hpXsc;
511 }
512 // pn to be np
513
514 if(neutron)
515 {
516 if( proj_momentum < 0.8 )
517 {
518 hnXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
519 }
520 else if( proj_momentum < 1.4 )
521 {
522 hnXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
523 }
524 else // if( proj_momentum < 10. )
525 {
526 hnXsc = 33.3+
527 20.8*(std::pow(proj_momentum,2.0)-1.35)/
528 (std::pow(proj_momentum,2.50)+0.95);
529 }
530 fTotalXsc = hnXsc;
531 }
532 // xsection = hpXsc*Zt + hnXsc*Nt;
533 // xsection = hpXsc*(Zt + Nt);
534 // xsection = hnXsc*(Zt + Nt);
535 }
536 // xsection *= 0.95;
537 }
538 else if(theParticle == theAProton && pORn)
539 {
540 if(proton)
541 {
542 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
543 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
544 }
545 if(proton)
546 {
547 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
548 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
549 }
550 fTotalXsc = xsection;
551 }
552 else if(theParticle == thePiPlus && pORn)
553 {
554 if(proton)
555 {
556 if(proj_momentum < 0.4)
557 {
558 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
559 hpXsc = Ex3+20.0;
560 }
561 else if(proj_momentum < 1.15)
562 {
563 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
564 hpXsc = Ex4+14.0;
565 }
566 else if(proj_momentum < 3.5)
567 {
568 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
569 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
570 hpXsc = Ex1+Ex2+27.5;
571 }
572 else // if(proj_momentum > 3.5) // mb
573 {
574 hpXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
575 }
576 fTotalXsc = hpXsc;
577 }
578
579// pi+n = pi-p??
580
581 if(neutron)
582 {
583 if(proj_momentum < 0.37)
584 {
585 hnXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
586 }
587 else if(proj_momentum<0.65)
588 {
589 hnXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
590 }
591 else if(proj_momentum<1.3)
592 {
593 hnXsc = 36.1+
594 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
595 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
596 }
597 else if(proj_momentum<3.0)
598 {
599 hnXsc = 36.1+0.079-4.313*std::log(proj_momentum)+
600 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
601 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
602 }
603 else // mb
604 {
605 hnXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
606 }
607 fTotalXsc = hnXsc;
608 }
609 // xsection = hpXsc*Zt + hnXsc*Nt;
610 }
611 else if(theParticle == thePiMinus && pORn)
612 {
613 // pi-n = pi+p??
614
615 if(neutron)
616 {
617 if(proj_momentum < 0.4)
618 {
619 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
620 hnXsc = Ex3+20.0;
621 }
622 else if(proj_momentum < 1.15)
623 {
624 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
625 hnXsc = Ex4+14.0;
626 }
627 else if(proj_momentum < 3.5)
628 {
629 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
630 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
631 hnXsc = Ex1+Ex2+27.5;
632 }
633 else // if(proj_momentum > 3.5) // mb
634 {
635 hnXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
636 }
637 fTotalXsc = hnXsc;
638 }
639 // pi-p
640
641 if(proton)
642 {
643 if(proj_momentum < 0.37)
644 {
645 hpXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
646 }
647 else if(proj_momentum<0.65)
648 {
649 hpXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
650 }
651 else if(proj_momentum<1.3)
652 {
653 hpXsc = 36.1+
654 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
655 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
656 }
657 else if(proj_momentum<3.0)
658 {
659 hpXsc = 36.1+0.079-4.313*std::log(proj_momentum)+
660 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
661 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
662 }
663 else // mb
664 {
665 hpXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
666 }
667 fTotalXsc = hpXsc;
668 }
669 // xsection = hpXsc*Zt + hnXsc*Nt;
670 }
671 else if(theParticle == theKPlus && pORn)
672 {
673 if(proton)
674 {
675 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
676 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
677 }
678 if(neutron)
679 {
680 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
681 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
682 }
683 fTotalXsc = xsection;
684 }
685 else if(theParticle == theKMinus && pORn)
686 {
687 if(proton)
688 {
689 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
690 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
691 }
692 if(neutron)
693 {
694 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
695 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
696 }
697 fTotalXsc = xsection;
698 }
699 else if(theParticle == theSMinus && pORn)
700 {
701 xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
702 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
703 }
704 else if(theParticle == theGamma && pORn) // modify later on
705 {
706 xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
707 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
708 fTotalXsc = xsection;
709 }
710 else // as proton ???
711 {
712 if(proton)
713 {
714 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
715 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
716 }
717 if(neutron)
718 {
719 xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
720 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
721 }
722 fTotalXsc = xsection;
723 }
724 fTotalXsc *= millibarn; // parametrised in mb
725 // xsection *= millibarn; // parametrised in mb
726
727 fInelasticXsc = 0.83*fTotalXsc;
728 fElasticXsc = fTotalXsc - fInelasticXsc;
729 if (fElasticXsc < 0.)fElasticXsc = 0.;
730
731 return fTotalXsc;
732}
733
734/////////////////////////////////////////////////////////////////////////////////////
735//
736// Returns hadron-nucleon cross-section based on V. Uzjinsky parametrisation of
737// data from G4FTFCrossSection class
738
739G4double
740G4HadronNucleonXsc::GetHadronNucleonXscVU(const G4DynamicParticle* aParticle,
741 const G4ParticleDefinition* nucleon )
742{
743 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
744 G4int absPDGcode = std::abs(PDGcode);
745 G4double Elab = aParticle->GetTotalEnergy();
746 // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
747 G4double Plab = aParticle->GetMomentum().mag();
748 // std::sqrt(Elab * Elab - 0.88);
749
750 Elab /= GeV;
751 Plab /= GeV;
752
753 G4double LogPlab = std::log( Plab );
754 G4double sqrLogPlab = LogPlab * LogPlab;
755
756 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
757 G4bool proton = (nucleon == theProton);
758 G4bool neutron = (nucleon == theNeutron);
759
760
761
762
763 if( absPDGcode > 1000 && pORn ) //------Projectile is baryon -
764 {
765 if(proton)
766 {
767 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
768 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
769 }
770 if(neutron)
771 {
772 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
773 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
774 }
775 }
776 else if( PDGcode == 211 && pORn ) //------Projectile is PionPlus ----
777 {
778 if(proton)
779 {
780 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
781 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
782 }
783 if(neutron)
784 {
785 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
786 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
787 }
788 }
789 else if( PDGcode == -211 && pORn ) //------Projectile is PionMinus ----
790 {
791 if(proton)
792 {
793 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
794 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
795 }
796 if(neutron)
797 {
798 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
799 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
800 }
801 }
802 else if( PDGcode == 111 && pORn ) //------Projectile is PionZero --
803 {
804 if(proton)
805 {
806 fTotalXsc = (16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
807 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
808
809 fElasticXsc = ( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
810 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
811
812 }
813 if(neutron)
814 {
815 fTotalXsc = (33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
816 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
817 fElasticXsc = ( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
818 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
819 }
820 }
821 else if( PDGcode == 321 && pORn ) //------Projectile is KaonPlus --
822 {
823 if(proton)
824 {
825 fTotalXsc = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
826 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
827 }
828 if(neutron)
829 {
830 fTotalXsc = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
831 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
832 }
833 }
834 else if( PDGcode ==-321 && pORn ) //------Projectile is KaonMinus ----
835 {
836 if(proton)
837 {
838 fTotalXsc = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66*sqrLogPlab - 5.6*LogPlab;
839 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29*sqrLogPlab - 2.4*LogPlab;
840 }
841 if(neutron)
842 {
843 fTotalXsc = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab;
844 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab;
845 }
846 }
847 else if( PDGcode == 311 && pORn ) //------Projectile is KaonZero -----
848 {
849 if(proton)
850 {
851 fTotalXsc = ( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
852 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
853 fElasticXsc = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
854 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
855 }
856 if(neutron)
857 {
858 fTotalXsc = ( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
859 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
860 fElasticXsc = ( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
861 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
862 }
863 }
864 else //------Projectile is undefined, Nucleon assumed
865 {
866 if(proton)
867 {
868 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
869 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
870 }
871 if(neutron)
872 {
873 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
874 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
875 }
876 }
877 fTotalXsc *= millibarn;
878 fElasticXsc *= millibarn;
879 fInelasticXsc = fTotalXsc - fElasticXsc;
880 if (fInelasticXsc < 0.) fInelasticXsc = 0.;
881
882 return fTotalXsc;
883}
884
885/////////////////////////////////////////////////////////////////////////////////////
886//
887// Returns hadron-nucleon cross-section based on Mikhail Kossov parametrisation of
888// data from G4QuasiFreeRatios class
889
890G4double
891G4HadronNucleonXsc::GetHadronNucleonXscMK(const G4DynamicParticle* aParticle,
892 const G4ParticleDefinition* nucleon )
893{
894 G4int I = -1;
895 G4int PDG = aParticle->GetDefinition()->GetPDGEncoding();
896
897 // G4int absPDG = std::abs(PDG);
898
899 G4double p = aParticle->GetMomentum().mag()/GeV;
900
901 G4bool F = false;
902 if(nucleon == theProton) F = true;
903 else if(nucleon == theNeutron) F = false;
904 else
905 {
906 G4cout << "nucleon is not proton or neutron, return xsc for proton" << G4endl;
907 F = true;
908 }
909
910 G4bool kfl = true; // Flag of K0/aK0 oscillation
911 G4bool kf = false;
912
913 if( PDG == 130 || PDG == 310 )
914 {
915 kf = true;
916 if( G4UniformRand() > .5 ) kfl = false;
917 }
918 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) I = 0; // pp/nn
919 else if( (PDG == 2112 && F) || (PDG == 2212 && !F) ) I = 1; // np/pn
920
921 else if( (PDG == -211 && F) || (PDG == 211 && !F) ) I = 2; // pimp/pipn
922 else if( (PDG == 211 && F) || (PDG ==-211 && !F) ) I = 3; // pipp/pimn
923
924 else if( PDG == -321 || PDG == -311 || ( kf && !kfl ) ) I = 4; // KmN/K0N
925 else if( PDG == 321 || PDG == 311 || ( kf && kfl ) ) I = 5; // KpN/aK0N
926
927 else if( PDG > 3000 && PDG < 3335) I = 6; // @@ for all hyperons - take Lambda
928 else if( PDG < -2000 && PDG > -3335) I = 7; // @@ for all anti-baryons - anti-p/anti-n
929 else
930 {
931 G4cout<<"MK PDG = "<<PDG
932 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
933 G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash");
934 }
935
936 // Each parameter set can have not more than nPoints = 128 parameters
937
938 static const G4double lmi = 3.5; // min of (lnP-lmi)^2 parabola
939 static const G4double pbe = .0557; // elastic (lnP-lmi)^2 parabola coefficient
940 static const G4double pbt = .3; // total (lnP-lmi)^2 parabola coefficient
941 static const G4double pmi = .1; // Below that fast LE calculation is made
942 static const G4double pma = 1000.; // Above that fast HE calculation is made
943
944 if( p <= 0.)
945 {
946 G4cout<<" p = "<<p<<" is zero or negative"<<G4endl;
947
948 fElasticXsc = 0.;
949 fInelasticXsc = 0.;
950 fTotalXsc = 0.;
951
952 return fTotalXsc;
953 }
954 if (!I) // pp/nn
955 {
956 if( p < pmi )
957 {
958 G4double p2 = p*p;
959 fElasticXsc = 1./(.00012 + p2*.2);
960 fTotalXsc = fElasticXsc;
961 }
962 else if(p>pma)
963 {
964 G4double lp = std::log(p)-lmi;
965 G4double lp2 = lp*lp;
966 fElasticXsc = pbe*lp2 + 6.72;
967 fTotalXsc = pbt*lp2 + 38.2;
968 }
969 else
970 {
971 G4double p2 = p*p;
972 G4double LE = 1./( .00012 + p2*.2);
973 G4double lp = std::log(p) - lmi;
974 G4double lp2 = lp*lp;
975 G4double rp2 = 1./p2;
976 fElasticXsc = LE + ( pbe*lp2 + 6.72+32.6/p)/( 1. + rp2/p);
977 fTotalXsc = LE + ( pbt*lp2 + 38.2+52.7*rp2)/( 1. + 2.72*rp2*rp2);
978 }
979 }
980 else if( I==1 ) // np/pn
981 {
982 if( p < pmi )
983 {
984 G4double p2 = p*p;
985 fElasticXsc = 1./( .00012 + p2*( .051 + .1*p2));
986 fTotalXsc = fElasticXsc;
987 }
988 else if( p > pma )
989 {
990 G4double lp = std::log(p) - lmi;
991 G4double lp2 = lp*lp;
992 fElasticXsc = pbe*lp2 + 6.72;
993 fTotalXsc = pbt*lp2 + 38.2;
994 }
995 else
996 {
997 G4double p2 = p*p;
998 G4double LE = 1./( .00012 + p2*( .051 + .1*p2 ) );
999 G4double lp = std::log(p) - lmi;
1000 G4double lp2 = lp*lp;
1001 G4double rp2 = 1./p2;
1002 fElasticXsc = LE + (pbe*lp2 + 6.72 + 30./p)/( 1. + .49*rp2/p);
1003 fTotalXsc = LE + (pbt*lp2 + 38.2)/( 1. + .54*rp2*rp2);
1004 }
1005 }
1006 else if( I == 2 ) // pimp/pipn
1007 {
1008 G4double lp = std::log(p);
1009
1010 if(p<pmi)
1011 {
1012 G4double lr = lp + 1.27;
1013 fElasticXsc = 1.53/( lr*lr + .0676);
1014 fTotalXsc = fElasticXsc*3;
1015 }
1016 else if( p > pma )
1017 {
1018 G4double ld = lp - lmi;
1019 G4double ld2 = ld*ld;
1020 G4double sp = std::sqrt(p);
1021 fElasticXsc = pbe*ld2 + 2.4 + 7./sp;
1022 fTotalXsc = pbt*ld2 + 22.3 + 12./sp;
1023 }
1024 else
1025 {
1026 G4double lr = lp + 1.27;
1027 G4double LE = 1.53/( lr*lr + .0676);
1028 G4double ld = lp - lmi;
1029 G4double ld2 = ld*ld;
1030 G4double p2 = p*p;
1031 G4double p4 = p2*p2;
1032 G4double sp = std::sqrt(p);
1033 G4double lm = lp + .36;
1034 G4double md = lm*lm + .04;
1035 G4double lh = lp - .017;
1036 G4double hd = lh*lh + .0025;
1037 fElasticXsc = LE + (pbe*ld2 + 2.4 + 7./sp)/( 1. + .7/p4) + .6/md + .05/hd;
1038 fTotalXsc = LE*3 + (pbt*ld2 + 22.3 + 12./sp)/(1. + .4/p4) + 1./md + .06/hd;
1039 }
1040 }
1041 else if( I == 3 ) // pipp/pimn
1042 {
1043 G4double lp = std::log(p);
1044
1045 if( p < pmi )
1046 {
1047 G4double lr = lp + 1.27;
1048 G4double lr2 = lr*lr;
1049 fElasticXsc = 13./( lr2 + lr2*lr2 + .0676);
1050 fTotalXsc = fElasticXsc;
1051 }
1052 else if( p > pma )
1053 {
1054 G4double ld = lp - lmi;
1055 G4double ld2 = ld*ld;
1056 G4double sp = std::sqrt(p);
1057 fElasticXsc = pbe*ld2 + 2.4 + 6./sp;
1058 fTotalXsc = pbt*ld2 + 22.3 + 5./sp;
1059 }
1060 else
1061 {
1062 G4double lr = lp + 1.27;
1063 G4double lr2 = lr*lr;
1064 G4double LE = 13./(lr2 + lr2*lr2 + .0676);
1065 G4double ld = lp - lmi;
1066 G4double ld2 = ld*ld;
1067 G4double p2 = p*p;
1068 G4double p4 = p2*p2;
1069 G4double sp = std::sqrt(p);
1070 G4double lm = lp - .32;
1071 G4double md = lm*lm + .0576;
1072 fElasticXsc = LE + (pbe*ld2 + 2.4 + 6./sp)/(1. + 3./p4) + .7/md;
1073 fTotalXsc = LE + (pbt*ld2 + 22.3 + 5./sp)/(1. + 1./p4) + .8/md;
1074 }
1075 }
1076 else if( I == 4 ) // Kmp/Kmn/K0p/K0n
1077 {
1078 if( p < pmi)
1079 {
1080 G4double psp = p*std::sqrt(p);
1081 fElasticXsc = 5.2/psp;
1082 fTotalXsc = 14./psp;
1083 }
1084 else if( p > pma )
1085 {
1086 G4double ld = std::log(p) - lmi;
1087 G4double ld2 = ld*ld;
1088 fElasticXsc = pbe*ld2 + 2.23;
1089 fTotalXsc = pbt*ld2 + 19.5;
1090 }
1091 else
1092 {
1093 G4double ld = std::log(p) - lmi;
1094 G4double ld2 = ld*ld;
1095 G4double sp = std::sqrt(p);
1096 G4double psp = p*sp;
1097 G4double p2 = p*p;
1098 G4double p4 = p2*p2;
1099 G4double lm = p - .39;
1100 G4double md = lm*lm + .000156;
1101 G4double lh = p - 1.;
1102 G4double hd = lh*lh + .0156;
1103 fElasticXsc = 5.2/psp + (pbe*ld2 + 2.23)/(1. - .7/sp + .075/p4) + .004/md + .15/hd;
1104 fTotalXsc = 14./psp + (pbt*ld2 + 19.5)/(1. - .21/sp + .52/p4) + .006/md + .30/hd;
1105 }
1106 }
1107 else if( I == 5 ) // Kpp/Kpn/aKp/aKn
1108 {
1109 if( p < pmi )
1110 {
1111 G4double lr = p - .38;
1112 G4double lm = p - 1.;
1113 G4double md = lm*lm + .372;
1114 fElasticXsc = .7/(lr*lr + .0676) + 2./md;
1115 fTotalXsc = fElasticXsc + .6/md;
1116 }
1117 else if( p > pma )
1118 {
1119 G4double ld = std::log(p) - lmi;
1120 G4double ld2 = ld*ld;
1121 fElasticXsc = pbe*ld2 + 2.23;
1122 fTotalXsc = pbt*ld2 + 19.5;
1123 }
1124 else
1125 {
1126 G4double ld = std::log(p) - lmi;
1127 G4double ld2 = ld*ld;
1128 G4double lr = p - .38;
1129 G4double LE = .7/(lr*lr + .0676);
1130 G4double sp = std::sqrt(p);
1131 G4double p2 = p*p;
1132 G4double p4 = p2*p2;
1133 G4double lm = p - 1.;
1134 G4double md = lm*lm + .372;
1135 fElasticXsc = LE + (pbe*ld2 + 2.23)/(1. - .7/sp + .1/p4) + 2./md;
1136 fTotalXsc = LE + (pbt*ld2 + 19.5)/(1. + .46/sp + 1.6/p4) + 2.6/md;
1137 }
1138 }
1139 else if( I == 6 ) // hyperon-N
1140 {
1141 if( p < pmi )
1142 {
1143 G4double p2 = p*p;
1144 fElasticXsc = 1./(.002 + p2*(.12 + p2));
1145 fTotalXsc = fElasticXsc;
1146 }
1147 else if( p > pma )
1148 {
1149 G4double lp = std::log(p) - lmi;
1150 G4double lp2 = lp*lp;
1151 G4double sp = std::sqrt(p);
1152 fElasticXsc = (pbe*lp2 + 6.72)/(1. + 2./sp);
1153 fTotalXsc = (pbt*lp2 + 38.2 + 900./sp)/(1. + 27./sp);
1154 }
1155 else
1156 {
1157 G4double p2 = p*p;
1158 G4double LE = 1./(.002 + p2*(.12 + p2));
1159 G4double lp = std::log(p) - lmi;
1160 G4double lp2 = lp*lp;
1161 G4double p4 = p2*p2;
1162 G4double sp = std::sqrt(p);
1163 fElasticXsc = LE + (pbe*lp2 + 6.72 + 99./p2)/(1. + 2./sp + 2./p4);
1164 fTotalXsc = LE + (pbt*lp2 + 38.2 + 900./sp)/(1. + 27./sp + 3./p4);
1165 }
1166 }
1167 else if( I == 7 ) // antibaryon-N
1168 {
1169 if( p > pma )
1170 {
1171 G4double lp = std::log(p) - lmi;
1172 G4double lp2 = lp*lp;
1173 fElasticXsc = pbe*lp2 + 6.72;
1174 fTotalXsc = pbt*lp2 + 38.2;
1175 }
1176 else
1177 {
1178 G4double ye = std::pow(p, 1.25);
1179 G4double yt = std::pow(p, .35);
1180 G4double lp = std::log(p) - lmi;
1181 G4double lp2 = lp*lp;
1182 fElasticXsc = 80./(ye + 1.) + pbe*lp2 + 6.72;
1183 fTotalXsc = (80./yt + .3)/yt +pbt*lp2 + 38.2;
1184 }
1185 }
1186 else
1187 {
1188 G4cout<<"PDG incoding = "<<I<<" is not defined (0-7)"<<G4endl;
1189
1190 }
1191 if( fElasticXsc > fTotalXsc ) fElasticXsc = fTotalXsc;
1192
1193 fTotalXsc *= millibarn;
1194 fElasticXsc *= millibarn;
1195 fInelasticXsc = fTotalXsc - fElasticXsc;
1196 if (fInelasticXsc < 0.) fInelasticXsc = 0.;
1197
1198
1199 return fTotalXsc;
1200}
1201
1202
1203
1204
1205
1206////////////////////////////////////////////////////////////////////////////////////
1207//
1208//
1209
1210G4double G4HadronNucleonXsc::CalculateEcmValue( const G4double mp ,
1211 const G4double mt ,
1212 const G4double Plab )
1213{
1214 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1215 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
1216 // G4double Pcm = Plab * mt / Ecm;
1217 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
1218
1219 return Ecm ; // KEcm;
1220}
1221
1222
1223////////////////////////////////////////////////////////////////////////////////////
1224//
1225//
1226
1227G4double G4HadronNucleonXsc::CalcMandelstamS( const G4double mp ,
1228 const G4double mt ,
1229 const G4double Plab )
1230{
1231 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
1232 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
1233
1234 return sMand;
1235}
1236
1237
1238//
1239//
1240///////////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.