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

Last change on this file since 1315 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 29.5 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//
888
889G4double G4HadronNucleonXsc::CalculateEcmValue( const G4double mp ,
890 const G4double mt ,
891 const G4double Plab )
892{
893 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
894 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
895 // G4double Pcm = Plab * mt / Ecm;
896 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
897
898 return Ecm ; // KEcm;
899}
900
901
902////////////////////////////////////////////////////////////////////////////////////
903//
904//
905
906G4double G4HadronNucleonXsc::CalcMandelstamS( const G4double mp ,
907 const G4double mt ,
908 const G4double Plab )
909{
910 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
911 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
912
913 return sMand;
914}
915
916
917//
918//
919///////////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.