source: trunk/source/processes/hadronic/cross_sections/src/G4GlauberGribovCrossSection.cc@ 830

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

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