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

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

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