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

Last change on this file was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 29.3 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// 14.03.07 V. Grichine - first implementation
27//
28
29#include "G4HadronNucleonXsc.hh"
30
31#include "G4ParticleTable.hh"
32#include "G4IonTable.hh"
33#include "G4ParticleDefinition.hh"
[1340]34#include "G4HadTmpUtil.hh"
[819]35
36
37G4HadronNucleonXsc::G4HadronNucleonXsc()
38: fUpperLimit( 10000 * GeV ),
[1347]39 fLowerLimit( 0.03 * MeV ),
40 fTotalXsc(0.0), fElasticXsc(0.0), fInelasticXsc(0.0), fHadronNucleonXsc(0.0)
[819]41{
42 theGamma = G4Gamma::Gamma();
43 theProton = G4Proton::Proton();
44 theNeutron = G4Neutron::Neutron();
45 theAProton = G4AntiProton::AntiProton();
46 theANeutron = G4AntiNeutron::AntiNeutron();
47 thePiPlus = G4PionPlus::PionPlus();
48 thePiMinus = G4PionMinus::PionMinus();
49 thePiZero = G4PionZero::PionZero();
50 theKPlus = G4KaonPlus::KaonPlus();
51 theKMinus = G4KaonMinus::KaonMinus();
52 theK0S = G4KaonZeroShort::KaonZeroShort();
53 theK0L = G4KaonZeroLong::KaonZeroLong();
54 theL = G4Lambda::Lambda();
55 theAntiL = G4AntiLambda::AntiLambda();
56 theSPlus = G4SigmaPlus::SigmaPlus();
57 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
58 theSMinus = G4SigmaMinus::SigmaMinus();
59 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
60 theS0 = G4SigmaZero::SigmaZero();
61 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
62 theXiMinus = G4XiMinus::XiMinus();
63 theXi0 = G4XiZero::XiZero();
64 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
65 theAXi0 = G4AntiXiZero::AntiXiZero();
66 theOmega = G4OmegaMinus::OmegaMinus();
67 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
68 theD = G4Deuteron::Deuteron();
69 theT = G4Triton::Triton();
70 theA = G4Alpha::Alpha();
71 theHe3 = G4He3::He3();
72}
73
74
75G4HadronNucleonXsc::~G4HadronNucleonXsc()
[1340]76{}
[819]77
78
79G4bool
80G4HadronNucleonXsc::IsApplicable(const G4DynamicParticle* aDP,
[1340]81 const G4Element* anElement)
[819]82{
[1340]83 G4int Z = G4lrint(anElement->GetZ());
84 G4int A = G4lrint(anElement->GetN());
85 return IsIsoApplicable(aDP, Z, A);
[819]86}
87
88////////////////////////////////////////////////////////////////////////////////////////
89//
90
91G4bool
[1340]92G4HadronNucleonXsc::IsIsoApplicable(const G4DynamicParticle* aDP,
93 G4int Z, G4int)
[819]94{
[1340]95 G4bool applicable = false;
[819]96 // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
97 G4double kineticEnergy = aDP->GetKineticEnergy();
98
99 const G4ParticleDefinition* theParticle = aDP->GetDefinition();
100
101 if ( ( kineticEnergy >= fLowerLimit &&
[1340]102 Z > 1 && // >= He
[819]103 ( theParticle == theAProton ||
104 theParticle == theGamma ||
105 theParticle == theKPlus ||
106 theParticle == theKMinus ||
107 theParticle == theSMinus) ) ||
108
109 ( kineticEnergy >= 0.1*fLowerLimit &&
[1340]110 Z > 1 && // >= He
[819]111 ( theParticle == theProton ||
112 theParticle == theNeutron ||
113 theParticle == thePiPlus ||
114 theParticle == thePiMinus ) ) ) applicable = true;
115
116 return applicable;
117}
118
119
120/////////////////////////////////////////////////////////////////////////////////////
121//
122// Returns hadron-nucleon Xsc according to differnt parametrisations:
123// [2] E. Levin, hep-ph/9710546
124// [3] U. Dersch, et al, hep-ex/9910052
125// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
126
127G4double
128G4HadronNucleonXsc::GetHadronNucleonXscEL(const G4DynamicParticle* aParticle,
129 const G4ParticleDefinition* nucleon )
130{
131 G4double xsection;
132
133
134 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
135
136 G4double proj_mass = aParticle->GetMass();
137 G4double proj_momentum = aParticle->GetMomentum().mag();
138 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
139
140 sMand /= GeV*GeV; // in GeV for parametrisation
141 proj_momentum /= GeV;
142
143 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
144
145 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
146
147
148 if(theParticle == theGamma && pORn )
149 {
150 xsection = (0.0677*std::pow(sMand,0.0808) + 0.129*std::pow(sMand,-0.4525));
151 }
152 else if(theParticle == theNeutron && pORn ) // as proton ???
153 {
154 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
155 }
156 else if(theParticle == theProton && pORn )
157 {
158 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
159
160 // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
161 // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
162 }
163 else if(theParticle == theAProton && pORn )
164 {
165 xsection = ( 21.70*std::pow(sMand,0.0808) + 98.39*std::pow(sMand,-0.4525));
166 }
167 else if(theParticle == thePiPlus && pORn )
168 {
169 xsection = (13.63*std::pow(sMand,0.0808) + 27.56*std::pow(sMand,-0.4525));
170 }
171 else if(theParticle == thePiMinus && pORn )
172 {
173 // xsection = At*( 55.2*std::pow(sMand,-0.255) + 0.346*std::log(sMand)*std::log(sMand) );
174 xsection = (13.63*std::pow(sMand,0.0808) + 36.02*std::pow(sMand,-0.4525));
175 }
176 else if(theParticle == theKPlus && pORn )
177 {
178 xsection = (11.82*std::pow(sMand,0.0808) + 8.15*std::pow(sMand,-0.4525));
179 }
180 else if(theParticle == theKMinus && pORn )
181 {
182 xsection = (11.82*std::pow(sMand,0.0808) + 26.36*std::pow(sMand,-0.4525));
183 }
184 else // as proton ???
185 {
186 xsection = (21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
187 }
188 xsection *= millibarn;
189
190 fTotalXsc = xsection;
191 fInelasticXsc = 0.83*xsection;
192 fElasticXsc = fTotalXsc - fInelasticXsc;
193 if (fElasticXsc < 0.)fElasticXsc = 0.;
194
195 return xsection;
196}
197
198
199/////////////////////////////////////////////////////////////////////////////////////
200//
201// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
202// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
203// At = number of nucleons, Zt = number of protons
204
205G4double
206G4HadronNucleonXsc::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
207 const G4ParticleDefinition* nucleon )
208{
209 G4double xsection(0);
210 G4int Zt=1, Nt=1, At=1;
211
212 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
213
214 G4double proj_mass = aParticle->GetMass();
215 G4double proj_momentum = aParticle->GetMomentum().mag();
216
217 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
218
219 sMand /= GeV*GeV; // in GeV for parametrisation
220
221 // General PDG fit constants
222
223 G4double s0 = 5.38*5.38; // in Gev^2
224 G4double eta1 = 0.458;
225 G4double eta2 = 0.458;
226 G4double B = 0.308;
227
228 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
229
230 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
231 G4bool proton = (nucleon == theProton);
232 G4bool neutron = (nucleon == theNeutron);
233
234 if(theParticle == theNeutron) // proton-neutron fit
235 {
236 if ( proton )
237 {
238 xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
239 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));// on p
240 }
241 if ( neutron )
242 {
243 xsection = Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
244 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // on n pp for nn
245 }
246 }
247 else if(theParticle == theProton)
248 {
249 if ( proton )
250 {
251 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
252 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
253 }
254 if ( neutron )
255 {
256 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
257 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
258 }
259 }
260 else if(theParticle == theAProton)
261 {
262 if ( proton )
263 {
264 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
265 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
266 }
267 if ( neutron )
268 {
269 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
270 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
271 }
272 }
273 else if(theParticle == thePiPlus && pORn )
274 {
275 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
276 + 19.24*std::pow(sMand,-eta1) - 6.03*std::pow(sMand,-eta2));
277 }
278 else if(theParticle == thePiMinus && pORn )
279 {
280 xsection = At*( 20.86 + B*std::pow(std::log(sMand/s0),2.)
281 + 19.24*std::pow(sMand,-eta1) + 6.03*std::pow(sMand,-eta2));
282 }
283 else if(theParticle == theKPlus)
284 {
285 if ( proton )
286 {
287 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
288 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
289 }
290 if ( neutron )
291 {
292 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
293 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
294 }
295 }
296 else if(theParticle == theKMinus)
297 {
298 if ( proton )
299 {
300 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
301 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
302 }
303 if ( neutron )
304 {
305 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
306 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2) );
307 }
308 }
309 else if(theParticle == theSMinus && pORn )
310 {
311 xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
312 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2) );
313 }
314 else if(theParticle == theGamma && pORn ) // modify later on
315 {
316 xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
317 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2) );
318
319 }
320 else // as proton ???
321 {
322 if ( proton )
323 {
324 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
325 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2) );
326 }
327 if ( neutron )
328 {
329 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
330 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
331 }
332 }
333 xsection *= millibarn; // parametrised in mb
334
335 fTotalXsc = xsection;
336 fInelasticXsc = 0.83*xsection;
337 fElasticXsc = fTotalXsc - fInelasticXsc;
338 if (fElasticXsc < 0.)fElasticXsc = 0.;
339
340 return xsection;
341}
342
343
344/////////////////////////////////////////////////////////////////////////////////////
345//
346// Returns hadron-nucleon cross-section based on N. Starkov parametrisation of
347// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
348
349G4double
350G4HadronNucleonXsc::GetHadronNucleonXscNS(const G4DynamicParticle* aParticle,
351 const G4ParticleDefinition* nucleon )
352{
353 G4double xsection(0), Delta, A0, B0;
354 G4int Zt=1, Nt=1, At=1;
355 G4double hpXsc(0);
356 G4double hnXsc(0);
357
358
359 G4double targ_mass = 0.939*GeV; // ~mean neutron and proton ???
360
361 G4double proj_mass = aParticle->GetMass();
362 G4double proj_energy = aParticle->GetTotalEnergy();
363 G4double proj_momentum = aParticle->GetMomentum().mag();
364
365 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
366
367 sMand /= GeV*GeV; // in GeV for parametrisation
368 proj_momentum /= GeV;
369 proj_energy /= GeV;
370 proj_mass /= GeV;
371
372 // General PDG fit constants
373
374 G4double s0 = 5.38*5.38; // in Gev^2
375 G4double eta1 = 0.458;
376 G4double eta2 = 0.458;
377 G4double B = 0.308;
378
379
380 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
381
382 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
383 G4bool proton = (nucleon == theProton);
384 G4bool neutron = (nucleon == theNeutron);
385
386 if( theParticle == theNeutron && pORn)
387 {
388 if( proj_momentum >= 10.)
389 // if( proj_momentum >= 2.)
390 {
391 Delta = 1.;
392
393 if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
394
395 if(proj_momentum >= 10.)
396 {
397 B0 = 7.5;
398 A0 = 100. - B0*std::log(3.0e7);
399
400 xsection = A0 + B0*std::log(proj_energy) - 11
401 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
402 0.93827*0.93827,-0.165); // mb
403 }
404 fTotalXsc = xsection;
405 }
406 else
407 {
408 // nn to be pp
409
410 if(neutron)
411 {
412 if( proj_momentum < 0.73 )
413 {
414 hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
415 }
416 else if( proj_momentum < 1.05 )
417 {
418 hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
419 (std::log(proj_momentum/0.73));
420 }
421 else // if( proj_momentum < 10. )
422 {
423 hnXsc = 39.0+
424 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
425 }
426 fTotalXsc = hnXsc;
427 }
428 // pn to be np
429
430 if(proton)
431 {
432 if( proj_momentum < 0.8 )
433 {
434 hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
435 }
436 else if( proj_momentum < 1.4 )
437 {
438 hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
439 }
440 else // if( proj_momentum < 10. )
441 {
442 hpXsc = 33.3+
443 20.8*(std::pow(proj_momentum,2.0)-1.35)/
444 (std::pow(proj_momentum,2.50)+0.95);
445 }
446 fTotalXsc = hpXsc;
447 }
448 // xsection = hpXsc*Zt + hnXsc*Nt;
449 }
450 }
451 else if(theParticle == theProton && pORn)
452 {
453 if( proj_momentum >= 10.)
454 // if( proj_momentum >= 2.)
455 {
456 Delta = 1.;
457
458 if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
459
460 if(proj_momentum >= 10.)
461 {
462 B0 = 7.5;
463 A0 = 100. - B0*std::log(3.0e7);
464
465 xsection = A0 + B0*std::log(proj_energy) - 11
466 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
467 0.93827*0.93827,-0.165); // mb
468 }
469 fTotalXsc = xsection;
470 }
471 else
472 {
473 // pp
474
475 if(proton)
476 {
477 if( proj_momentum < 0.73 )
478 {
479 hpXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
480 }
481 else if( proj_momentum < 1.05 )
482 {
483 hpXsc = 23 + 40*(std::log(proj_momentum/0.73))*
484 (std::log(proj_momentum/0.73));
485 }
486 else // if( proj_momentum < 10. )
487 {
488 hpXsc = 39.0+
489 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
490 }
491 fTotalXsc = hpXsc;
492 }
493 // pn to be np
494
495 if(neutron)
496 {
497 if( proj_momentum < 0.8 )
498 {
499 hnXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
500 }
501 else if( proj_momentum < 1.4 )
502 {
503 hnXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
504 }
505 else // if( proj_momentum < 10. )
506 {
507 hnXsc = 33.3+
508 20.8*(std::pow(proj_momentum,2.0)-1.35)/
509 (std::pow(proj_momentum,2.50)+0.95);
510 }
511 fTotalXsc = hnXsc;
512 }
513 // xsection = hpXsc*Zt + hnXsc*Nt;
514 // xsection = hpXsc*(Zt + Nt);
515 // xsection = hnXsc*(Zt + Nt);
516 }
517 // xsection *= 0.95;
518 }
519 else if(theParticle == theAProton && pORn)
520 {
521 if(proton)
522 {
523 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
524 + 42.53*std::pow(sMand,-eta1) + 33.34*std::pow(sMand,-eta2));
525 }
526 if(proton)
527 {
528 xsection = Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
529 + 40.15*std::pow(sMand,-eta1) + 30.*std::pow(sMand,-eta2));
530 }
531 fTotalXsc = xsection;
532 }
533 else if(theParticle == thePiPlus && pORn)
534 {
535 if(proton)
536 {
537 if(proj_momentum < 0.4)
538 {
539 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
540 hpXsc = Ex3+20.0;
541 }
542 else if(proj_momentum < 1.15)
543 {
544 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
545 hpXsc = Ex4+14.0;
546 }
547 else if(proj_momentum < 3.5)
548 {
549 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
550 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
551 hpXsc = Ex1+Ex2+27.5;
552 }
553 else // if(proj_momentum > 3.5) // mb
554 {
555 hpXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
556 }
557 fTotalXsc = hpXsc;
558 }
559
560// pi+n = pi-p??
561
562 if(neutron)
563 {
564 if(proj_momentum < 0.37)
565 {
566 hnXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
567 }
568 else if(proj_momentum<0.65)
569 {
570 hnXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
571 }
572 else if(proj_momentum<1.3)
573 {
574 hnXsc = 36.1+
575 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
576 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
577 }
578 else if(proj_momentum<3.0)
579 {
580 hnXsc = 36.1+0.079-4.313*std::log(proj_momentum)+
581 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
582 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
583 }
584 else // mb
585 {
586 hnXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
587 }
588 fTotalXsc = hnXsc;
589 }
590 // xsection = hpXsc*Zt + hnXsc*Nt;
591 }
592 else if(theParticle == thePiMinus && pORn)
593 {
594 // pi-n = pi+p??
595
596 if(neutron)
597 {
598 if(proj_momentum < 0.4)
599 {
600 G4double Ex3 = 180*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.085/0.085);
601 hnXsc = Ex3+20.0;
602 }
603 else if(proj_momentum < 1.15)
604 {
605 G4double Ex4 = 88*(std::log(proj_momentum/0.75))*(std::log(proj_momentum/0.75));
606 hnXsc = Ex4+14.0;
607 }
608 else if(proj_momentum < 3.5)
609 {
610 G4double Ex1 = 3.2*std::exp(-(proj_momentum-2.55)*(proj_momentum-2.55)/0.55/0.55);
611 G4double Ex2 = 12*std::exp(-(proj_momentum-1.47)*(proj_momentum-1.47)/0.225/0.225);
612 hnXsc = Ex1+Ex2+27.5;
613 }
614 else // if(proj_momentum > 3.5) // mb
615 {
616 hnXsc = 10.6+2.*std::log(proj_energy)+25*std::pow(proj_energy,-0.43);
617 }
618 fTotalXsc = hnXsc;
619 }
620 // pi-p
621
622 if(proton)
623 {
624 if(proj_momentum < 0.37)
625 {
626 hpXsc = 28.0 + 40*std::exp(-(proj_momentum-0.29)*(proj_momentum-0.29)/0.07/0.07);
627 }
628 else if(proj_momentum<0.65)
629 {
630 hpXsc = 26+110*(std::log(proj_momentum/0.48))*(std::log(proj_momentum/0.48));
631 }
632 else if(proj_momentum<1.3)
633 {
634 hpXsc = 36.1+
635 10*std::exp(-(proj_momentum-0.72)*(proj_momentum-0.72)/0.06/0.06)+
636 24*std::exp(-(proj_momentum-1.015)*(proj_momentum-1.015)/0.075/0.075);
637 }
638 else if(proj_momentum<3.0)
639 {
640 hpXsc = 36.1+0.079-4.313*std::log(proj_momentum)+
641 3*std::exp(-(proj_momentum-2.1)*(proj_momentum-2.1)/0.4/0.4)+
642 1.5*std::exp(-(proj_momentum-1.4)*(proj_momentum-1.4)/0.12/0.12);
643 }
644 else // mb
645 {
646 hpXsc = 10.6+2*std::log(proj_energy)+30*std::pow(proj_energy,-0.43);
647 }
648 fTotalXsc = hpXsc;
649 }
650 // xsection = hpXsc*Zt + hnXsc*Nt;
651 }
652 else if(theParticle == theKPlus && pORn)
653 {
654 if(proton)
655 {
656 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
657 + 7.14*std::pow(sMand,-eta1) - 13.45*std::pow(sMand,-eta2));
658 }
659 if(neutron)
660 {
661 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
662 + 5.17*std::pow(sMand,-eta1) - 7.23*std::pow(sMand,-eta2));
663 }
664 fTotalXsc = xsection;
665 }
666 else if(theParticle == theKMinus && pORn)
667 {
668 if(proton)
669 {
670 xsection = Zt*( 17.91 + B*std::pow(std::log(sMand/s0),2.)
671 + 7.14*std::pow(sMand,-eta1) + 13.45*std::pow(sMand,-eta2));
672 }
673 if(neutron)
674 {
675 xsection = Nt*( 17.87 + B*std::pow(std::log(sMand/s0),2.)
676 + 5.17*std::pow(sMand,-eta1) + 7.23*std::pow(sMand,-eta2));
677 }
678 fTotalXsc = xsection;
679 }
680 else if(theParticle == theSMinus && pORn)
681 {
682 xsection = At*( 35.20 + B*std::pow(std::log(sMand/s0),2.)
683 - 199.*std::pow(sMand,-eta1) + 264.*std::pow(sMand,-eta2));
684 }
685 else if(theParticle == theGamma && pORn) // modify later on
686 {
687 xsection = At*( 0.0 + B*std::pow(std::log(sMand/s0),2.)
688 + 0.032*std::pow(sMand,-eta1) - 0.0*std::pow(sMand,-eta2));
689 fTotalXsc = xsection;
690 }
691 else // as proton ???
692 {
693 if(proton)
694 {
695 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
696 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
697 }
698 if(neutron)
699 {
700 xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
701 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
702 }
703 fTotalXsc = xsection;
704 }
705 fTotalXsc *= millibarn; // parametrised in mb
706 // xsection *= millibarn; // parametrised in mb
707
708 fInelasticXsc = 0.83*fTotalXsc;
709 fElasticXsc = fTotalXsc - fInelasticXsc;
710 if (fElasticXsc < 0.)fElasticXsc = 0.;
711
712 return fTotalXsc;
713}
714
715/////////////////////////////////////////////////////////////////////////////////////
716//
717// Returns hadron-nucleon cross-section based on V. Uzjinsky parametrisation of
718// data from G4FTFCrossSection class
719
720G4double
721G4HadronNucleonXsc::GetHadronNucleonXscVU(const G4DynamicParticle* aParticle,
722 const G4ParticleDefinition* nucleon )
723{
724 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
725 G4int absPDGcode = std::abs(PDGcode);
726 G4double Elab = aParticle->GetTotalEnergy();
727 // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
728 G4double Plab = aParticle->GetMomentum().mag();
729 // std::sqrt(Elab * Elab - 0.88);
730
731 Elab /= GeV;
732 Plab /= GeV;
733
734 G4double LogPlab = std::log( Plab );
735 G4double sqrLogPlab = LogPlab * LogPlab;
736
737 G4bool pORn = (nucleon == theProton || nucleon == theNeutron );
738 G4bool proton = (nucleon == theProton);
739 G4bool neutron = (nucleon == theNeutron);
740
741
742 if( absPDGcode > 1000 && pORn ) //------Projectile is baryon -
743 {
744 if(proton)
745 {
746 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
747 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
748 }
749 if(neutron)
750 {
751 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
752 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
753 }
754 }
755 else if( PDGcode == 211 && pORn ) //------Projectile is PionPlus ----
756 {
757 if(proton)
758 {
759 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
760 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
761 }
762 if(neutron)
763 {
764 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
765 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
766 }
767 }
768 else if( PDGcode == -211 && pORn ) //------Projectile is PionMinus ----
769 {
770 if(proton)
771 {
772 fTotalXsc = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
773 fElasticXsc = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
774 }
775 if(neutron)
776 {
777 fTotalXsc = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
778 fElasticXsc = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
779 }
780 }
781 else if( PDGcode == 111 && pORn ) //------Projectile is PionZero --
782 {
783 if(proton)
784 {
785 fTotalXsc = (16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
786 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
787
788 fElasticXsc = ( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
789 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
790
791 }
792 if(neutron)
793 {
794 fTotalXsc = (33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
795 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
796 fElasticXsc = ( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
797 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
798 }
799 }
800 else if( PDGcode == 321 && pORn ) //------Projectile is KaonPlus --
801 {
802 if(proton)
803 {
804 fTotalXsc = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
805 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
806 }
807 if(neutron)
808 {
809 fTotalXsc = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
810 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
811 }
812 }
813 else if( PDGcode ==-321 && pORn ) //------Projectile is KaonMinus ----
814 {
815 if(proton)
816 {
817 fTotalXsc = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66*sqrLogPlab - 5.6*LogPlab;
818 fElasticXsc = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29*sqrLogPlab - 2.4*LogPlab;
819 }
820 if(neutron)
821 {
822 fTotalXsc = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38*sqrLogPlab - 2.9*LogPlab;
823 fElasticXsc = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16*sqrLogPlab - 1.3*LogPlab;
824 }
825 }
826 else if( PDGcode == 311 && pORn ) //------Projectile is KaonZero -----
827 {
828 if(proton)
829 {
830 fTotalXsc = ( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
831 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
832 fElasticXsc = ( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
833 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
834 }
835 if(neutron)
836 {
837 fTotalXsc = ( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
838 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
839 fElasticXsc = ( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
840 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
841 }
842 }
843 else //------Projectile is undefined, Nucleon assumed
844 {
845 if(proton)
846 {
847 fTotalXsc = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
848 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
849 }
850 if(neutron)
851 {
852 fTotalXsc = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
853 fElasticXsc = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
854 }
855 }
856 fTotalXsc *= millibarn;
857 fElasticXsc *= millibarn;
858 fInelasticXsc = fTotalXsc - fElasticXsc;
859 if (fInelasticXsc < 0.) fInelasticXsc = 0.;
860
861 return fTotalXsc;
862}
863
864////////////////////////////////////////////////////////////////////////////////////
865//
866//
867
868G4double G4HadronNucleonXsc::CalculateEcmValue( const G4double mp ,
869 const G4double mt ,
870 const G4double Plab )
871{
872 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
873 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
874 // G4double Pcm = Plab * mt / Ecm;
875 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
876
877 return Ecm ; // KEcm;
878}
879
880
881////////////////////////////////////////////////////////////////////////////////////
882//
883//
884
885G4double G4HadronNucleonXsc::CalcMandelstamS( const G4double mp ,
886 const G4double mt ,
887 const G4double Plab )
888{
889 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
890 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
891
892 return sMand;
893}
894
895
896//
897//
898///////////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.