source: trunk/source/processes/hadronic/cross_sections/src/G4GGNuclNuclCrossSection.cc@ 1036

Last change on this file since 1036 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 26.6 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// 24.11.08 V. Grichine - first implementation
28//
29
30#include "G4GGNuclNuclCrossSection.hh"
31
32#include "G4ParticleTable.hh"
33#include "G4IonTable.hh"
34#include "G4ParticleDefinition.hh"
35
36
37
38////////////////////////////////////////////////////////////////////////////////
39//
40//
41
42G4GGNuclNuclCrossSection::G4GGNuclNuclCrossSection()
43: fUpperLimit( 100000 * GeV ),
44 fLowerLimit( 0.1 * GeV ),
45 fRadiusConst( 1.08*fermi ) // 1.1, 1.3 ?
46{
47 theProton = G4Proton::Proton();
48 theNeutron = G4Neutron::Neutron();
49}
50
51///////////////////////////////////////////////////////////////////////////////////////
52//
53//
54
55G4GGNuclNuclCrossSection::~G4GGNuclNuclCrossSection()
56{
57}
58
59
60////////////////////////////////////////////////////////////////////////////////////////
61//
62//
63
64
65G4bool
66G4GGNuclNuclCrossSection::IsApplicable(const G4DynamicParticle* aDP,
67 const G4Element* anElement)
68{
69 return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN());
70}
71
72////////////////////////////////////////////////////////////////////////////////////////
73//
74//
75
76G4bool
77G4GGNuclNuclCrossSection::IsZAApplicable(const G4DynamicParticle* aDP,
78 G4double Z, G4double)
79{
80 G4bool applicable = false;
81 // G4int baryonNumber = aDP->GetDefinition()->GetBaryonNumber();
82 G4double kineticEnergy = aDP->GetKineticEnergy();
83
84 // const G4ParticleDefinition* theParticle = aDP->GetDefinition();
85
86 if ( kineticEnergy >= fLowerLimit && Z > 1.5 ) applicable = true;
87
88 return applicable;
89}
90
91////////////////////////////////////////////////////////////////////////////////////////
92//
93// Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
94// Glauber model with Gribov correction calculated in the dipole approximation on
95// light cone. Gaussian density helps to calculate rest integrals of the model.
96// [1] B.Z. Kopeliovich, nucl-th/0306044
97
98
99G4double G4GGNuclNuclCrossSection::
100GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement, G4double T)
101{
102 return GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(), T);
103}
104
105////////////////////////////////////////////////////////////////////////////////////////
106//
107// Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
108// Glauber model with Gribov correction calculated in the dipole approximation on
109// light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model.
110// [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
111
112
113
114G4double G4GGNuclNuclCrossSection::
115GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double tZ, G4double tA, G4double)
116{
117 G4double xsection, sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
118
119 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
120 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
121
122 G4double pTkin = aParticle->GetKineticEnergy();
123 pTkin /= pA;
124
125 G4double pN = pA - pZ;
126 if( pN < 0. ) pN = 0.;
127
128 G4double tN = tA - tZ;
129 if( tN < 0. ) tN = 0.;
130
131 G4double tR = GetNucleusRadius(tA);
132 G4double pR = GetNucleusRadius(pA);
133
134 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
135 (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
136
137 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
138
139 ratio = sigma/nucleusSquare;
140
141 xsection = nucleusSquare*std::log( 1. + ratio );
142
143 fTotalXsc = xsection;
144
145 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
146
147 fElasticXsc = fTotalXsc - fInelasticXsc;
148
149 G4double difratio = ratio/(1.+ratio);
150
151 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
152
153 // production to be checked !!! edit MK xsc
154
155 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscMK(theProton, pTkin, theProton) +
156 (pZ*tN+pN*tZ)*GetHadronNucleonXscMK(theProton, pTkin, theNeutron);
157
158 ratio = sigma/nucleusSquare;
159
160 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
161
162 if (fElasticXsc < 0.) fElasticXsc = 0.;
163
164 return xsection;
165}
166
167//////////////////////////////////////////////////////////////////////////
168//
169// Return single-diffraction/inelastic cross-section ratio
170
171G4double G4GGNuclNuclCrossSection::
172GetRatioSD(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
173{
174 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
175
176 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
177 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
178
179 G4double pTkin = aParticle->GetKineticEnergy();
180 pTkin /= pA;
181
182 G4double pN = pA - pZ;
183 if( pN < 0. ) pN = 0.;
184
185 G4double tN = tA - tZ;
186 if( tN < 0. ) tN = 0.;
187
188 G4double tR = GetNucleusRadius(tA);
189 G4double pR = GetNucleusRadius(pA);
190
191 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
192 (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
193
194 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
195
196 ratio = sigma/nucleusSquare;
197
198
199 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
200
201 G4double difratio = ratio/(1.+ratio);
202
203 fDiffractionXsc = 0.5*nucleusSquare*( difratio - std::log( 1. + difratio ) );
204
205 if (fInelasticXsc > 0.) ratio = fDiffractionXsc/fInelasticXsc;
206 else ratio = 0.;
207
208 return ratio;
209}
210
211//////////////////////////////////////////////////////////////////////////
212//
213// Return suasi-elastic/inelastic cross-section ratio
214
215G4double G4GGNuclNuclCrossSection::
216GetRatioQE(const G4DynamicParticle* aParticle, G4double tA, G4double tZ)
217{
218 G4double sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, ratio;
219
220 G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
221 G4double pA = aParticle->GetDefinition()->GetBaryonNumber();
222
223 G4double pTkin = aParticle->GetKineticEnergy();
224 pTkin /= pA;
225
226 G4double pN = pA - pZ;
227 if( pN < 0. ) pN = 0.;
228
229 G4double tN = tA - tZ;
230 if( tN < 0. ) tN = 0.;
231
232 G4double tR = GetNucleusRadius(tA);
233 G4double pR = GetNucleusRadius(pA);
234
235 sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
236 (pZ*tN+pN*tZ)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
237
238 nucleusSquare = cofTotal*pi*( pR*pR + tR*tR ); // basically 2piRR
239
240 ratio = sigma/nucleusSquare;
241
242 fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
243
244 // sigma = GetHNinelasticXsc(aParticle, tA, tZ);
245 ratio = sigma/nucleusSquare;
246
247 fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
248
249 if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
250 else ratio = 0.;
251 if ( ratio < 0. ) ratio = 0.;
252
253 return ratio;
254}
255
256/////////////////////////////////////////////////////////////////////////////////////
257//
258// Returns hadron-nucleon Xsc according to differnt parametrisations:
259// [2] E. Levin, hep-ph/9710546
260// [3] U. Dersch, et al, hep-ex/9910052
261// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
262
263G4double
264G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
265 const G4Element* anElement )
266{
267 G4double At = anElement->GetN(); // number of nucleons
268 G4double Zt = anElement->GetZ(); // number of protons
269
270
271 return GetHadronNucleonXsc( aParticle, At, Zt );
272}
273
274
275
276
277/////////////////////////////////////////////////////////////////////////////////////
278//
279// Returns hadron-nucleon Xsc according to differnt parametrisations:
280// [2] E. Levin, hep-ph/9710546
281// [3] U. Dersch, et al, hep-ex/9910052
282// [4] M.J. Longo, et al, Phys.Rev.Lett. 33 (1974) 725
283
284G4double
285G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
286 G4double At, G4double Zt )
287{
288 G4double xsection = 0.;
289
290
291 G4double targ_mass = G4ParticleTable::GetParticleTable()->
292 GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
293
294 targ_mass = 0.939*GeV; // ~mean neutron and proton ???
295
296 G4double proj_mass = aParticle->GetMass();
297 G4double proj_momentum = aParticle->GetMomentum().mag();
298 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
299
300 sMand /= GeV*GeV; // in GeV for parametrisation
301 proj_momentum /= GeV;
302
303 const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
304
305
306 if(pParticle == theNeutron) // as proton ???
307 {
308 xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
309 }
310 else if(pParticle == theProton)
311 {
312 xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
313 // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
314 // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
315 }
316
317 xsection *= millibarn;
318 return xsection;
319}
320
321
322/////////////////////////////////////////////////////////////////////////////////////
323//
324// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
325// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
326
327G4double
328G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
329 const G4Element* anElement )
330{
331 G4double At = anElement->GetN(); // number of nucleons
332 G4double Zt = anElement->GetZ(); // number of protons
333
334
335 return GetHadronNucleonXscPDG( aParticle, At, Zt );
336}
337
338
339
340
341/////////////////////////////////////////////////////////////////////////////////////
342//
343// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
344// http://pdg.lbl.gov/2006/reviews/hadronicrpp.pdf
345// At = number of nucleons, Zt = number of protons
346
347G4double
348G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
349 G4double At, G4double Zt )
350{
351 G4double xsection = 0.;
352
353 G4double Nt = At-Zt; // number of neutrons
354 if (Nt < 0.) Nt = 0.;
355
356
357 G4double targ_mass = G4ParticleTable::GetParticleTable()->
358 GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
359
360 targ_mass = 0.939*GeV; // ~mean neutron and proton ???
361
362 G4double proj_mass = aParticle->GetMass();
363 G4double proj_momentum = aParticle->GetMomentum().mag();
364
365 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
366
367 sMand /= GeV*GeV; // in GeV for parametrisation
368
369 // General PDG fit constants
370
371 G4double s0 = 5.38*5.38; // in Gev^2
372 G4double eta1 = 0.458;
373 G4double eta2 = 0.458;
374 G4double B = 0.308;
375
376
377 const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
378
379
380 if(pParticle == theNeutron) // proton-neutron fit
381 {
382 xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
383 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
384 xsection += Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
385 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
386 }
387 else if(pParticle == theProton)
388 {
389
390 xsection = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
391 + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
392
393 xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
394 + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
395 }
396 xsection *= millibarn; // parametrised in mb
397 return xsection;
398}
399
400
401
402
403/////////////////////////////////////////////////////////////////////////////////////
404//
405// Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
406// data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
407// projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
408
409G4double
410G4GGNuclNuclCrossSection::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
411 G4double pTkin,
412 G4ParticleDefinition* tParticle)
413{
414 G4double xsection(0), Delta, A0, B0;
415 G4double hpXsc(0);
416 G4double hnXsc(0);
417
418
419 G4double targ_mass = tParticle->GetPDGMass();
420 G4double proj_mass = pParticle->GetPDGMass();
421
422 G4double proj_energy = proj_mass + pTkin;
423 G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
424
425 G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
426
427 sMand /= GeV*GeV; // in GeV for parametrisation
428 proj_momentum /= GeV;
429 proj_energy /= GeV;
430 proj_mass /= GeV;
431
432 // General PDG fit constants
433
434 // G4double s0 = 5.38*5.38; // in Gev^2
435 // G4double eta1 = 0.458;
436 // G4double eta2 = 0.458;
437 // G4double B = 0.308;
438
439
440
441
442
443 if( proj_momentum >= 10. ) // high energy: pp = nn = np
444 // if( proj_momentum >= 2.)
445 {
446 Delta = 1.;
447
448 if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
449
450 if( proj_momentum >= 10.)
451 {
452 B0 = 7.5;
453 A0 = 100. - B0*std::log(3.0e7);
454
455 xsection = A0 + B0*std::log(proj_energy) - 11
456 + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
457 0.93827*0.93827,-0.165); // mb
458 }
459 }
460 else // low energy pp = nn != np
461 {
462 if(pParticle == tParticle) // pp or nn // nn to be pp
463 {
464 if( proj_momentum < 0.73 )
465 {
466 hnXsc = 23 + 50*( std::pow( std::log(0.73/proj_momentum), 3.5 ) );
467 }
468 else if( proj_momentum < 1.05 )
469 {
470 hnXsc = 23 + 40*(std::log(proj_momentum/0.73))*
471 (std::log(proj_momentum/0.73));
472 }
473 else // if( proj_momentum < 10. )
474 {
475 hnXsc = 39.0 +
476 75*(proj_momentum - 1.2)/(std::pow(proj_momentum,3.0) + 0.15);
477 }
478 xsection = hnXsc;
479 }
480 else // pn to be np
481 {
482 if( proj_momentum < 0.8 )
483 {
484 hpXsc = 33+30*std::pow(std::log(proj_momentum/1.3),4.0);
485 }
486 else if( proj_momentum < 1.4 )
487 {
488 hpXsc = 33+30*std::pow(std::log(proj_momentum/0.95),2.0);
489 }
490 else // if( proj_momentum < 10. )
491 {
492 hpXsc = 33.3+
493 20.8*(std::pow(proj_momentum,2.0)-1.35)/
494 (std::pow(proj_momentum,2.50)+0.95);
495 }
496 xsection = hpXsc;
497 }
498 }
499 xsection *= millibarn; // parametrised in mb
500 return xsection;
501}
502
503/*
504/////////////////////////////////////////////////////////////////////////////////////
505//
506// Returns hadron-nucleon inelastic cross-section based on proper parametrisation
507
508G4double
509G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
510 const G4Element* anElement )
511{
512 G4double At = anElement->GetN(); // number of nucleons
513 G4double Zt = anElement->GetZ(); // number of protons
514
515
516 return GetHNinelasticXsc( aParticle, At, Zt );
517}
518
519/////////////////////////////////////////////////////////////////////////////////////
520//
521// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
522
523G4double
524G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
525 G4double At, G4double Zt )
526{
527 // G4ParticleDefinition* hadron = aParticle->GetDefinition();
528 G4double sumInelastic, Nt = At - Zt;
529
530 if(Nt < 0.) Nt = 0.;
531
532 sumInelastic = Zt*GetHadronNucleonXscMK(aParticle, theProton);
533 sumInelastic += Nt*GetHadronNucleonXscMK(aParticle, theNeutron);
534
535 return sumInelastic;
536}
537*/
538
539/////////////////////////////////////////////////////////////////////////////////////
540//
541// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
542
543G4double
544G4GGNuclNuclCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
545 G4double At, G4double Zt )
546{
547 G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
548 G4int absPDGcode = std::abs(PDGcode);
549
550 G4double Elab = aParticle->GetTotalEnergy();
551 // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
552 G4double Plab = aParticle->GetMomentum().mag();
553 // std::sqrt(Elab * Elab - 0.88);
554
555 Elab /= GeV;
556 Plab /= GeV;
557
558 G4double LogPlab = std::log( Plab );
559 G4double sqrLogPlab = LogPlab * LogPlab;
560
561 //G4cout<<"Plab = "<<Plab<<G4endl;
562
563 G4double NumberOfTargetProtons = Zt;
564 G4double NumberOfTargetNucleons = At;
565 G4double NumberOfTargetNeutrons = NumberOfTargetNucleons - NumberOfTargetProtons;
566
567 if(NumberOfTargetNeutrons < 0.) NumberOfTargetNeutrons = 0.;
568
569 G4double Xtotal = 0., Xelastic = 0., Xinelastic =0.;
570
571 if( absPDGcode > 1000 ) //------Projectile is baryon --------
572 {
573 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) +
574 0.522*sqrLogPlab - 4.51*LogPlab;
575
576 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) +
577 0.513*sqrLogPlab - 4.27*LogPlab;
578
579 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) +
580 0.169*sqrLogPlab - 1.85*LogPlab;
581
582 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) +
583 0.169*sqrLogPlab - 1.85*LogPlab;
584
585 Xtotal = ( NumberOfTargetProtons * XtotPP +
586 NumberOfTargetNeutrons * XtotPN );
587
588 Xelastic = ( NumberOfTargetProtons * XelPP +
589 NumberOfTargetNeutrons * XelPN );
590 }
591 Xinelastic = Xtotal - Xelastic;
592
593 if(Xinelastic < 0.) Xinelastic = 0.;
594
595 return Xinelastic*= millibarn;
596}
597
598/////////////////////////////////////////////////////////////////////////////////////
599//
600// Returns hadron-nucleon cross-section based on Mikhail Kossov CHIPS parametrisation of
601// data from G4QuasiFreeRatios class
602
603G4double
604G4GGNuclNuclCrossSection::GetHadronNucleonXscMK(G4ParticleDefinition* pParticle, G4double pTkin,
605 G4ParticleDefinition* nucleon )
606{
607 G4int I = -1;
608 G4int PDG = pParticle->GetPDGEncoding();
609 G4double totalXsc = 0;
610 G4double elasticXsc = 0;
611 G4double inelasticXsc;
612 // G4int absPDG = std::abs(PDG);
613
614 G4double pM = pParticle->GetPDGMass();
615 G4double p = std::sqrt(pTkin*(pTkin+2*pM))/GeV;
616
617 G4bool F = false;
618 if(nucleon == theProton) F = true;
619 else if(nucleon == theNeutron) F = false;
620 else
621 {
622 G4cout << "nucleon is not proton or neutron, return xsc for proton" << G4endl;
623 F = true;
624 }
625
626 G4bool kfl = true; // Flag of K0/aK0 oscillation
627 G4bool kf = false;
628
629 if( PDG == 130 || PDG == 310 )
630 {
631 kf = true;
632 if( G4UniformRand() > .5 ) kfl = false;
633 }
634 if ( (PDG == 2212 && F) || (PDG == 2112 && !F) ) I = 0; // pp/nn
635 else if( (PDG == 2112 && F) || (PDG == 2212 && !F) ) I = 1; // np/pn
636 else
637 {
638 G4cout<<"MK PDG = "<<PDG
639 <<", while it is defined only for p,n,hyperons,anti-baryons,pi,K/antiK"<<G4endl;
640 G4Exception("G4QuasiFreeRatio::FetchElTot:","22",FatalException,"CHIPScrash");
641 }
642
643 // Each parameter set can have not more than nPoints = 128 parameters
644
645 static const G4double lmi = 3.5; // min of (lnP-lmi)^2 parabola
646 static const G4double pbe = .0557; // elastic (lnP-lmi)^2 parabola coefficient
647 static const G4double pbt = .3; // total (lnP-lmi)^2 parabola coefficient
648 static const G4double pmi = .1; // Below that fast LE calculation is made
649 static const G4double pma = 1000.; // Above that fast HE calculation is made
650
651 if( p <= 0.)
652 {
653 G4cout<<" p = "<<p<<" is zero or negative"<<G4endl;
654
655 elasticXsc = 0.;
656 inelasticXsc = 0.;
657 totalXsc = 0.;
658
659 return totalXsc;
660 }
661 if (!I) // pp/nn
662 {
663 if( p < pmi )
664 {
665 G4double p2 = p*p;
666 elasticXsc = 1./(.00012 + p2*.2);
667 totalXsc = elasticXsc;
668 }
669 else if(p>pma)
670 {
671 G4double lp = std::log(p)-lmi;
672 G4double lp2 = lp*lp;
673 elasticXsc = pbe*lp2 + 6.72;
674 totalXsc = pbt*lp2 + 38.2;
675 }
676 else
677 {
678 G4double p2 = p*p;
679 G4double LE = 1./( .00012 + p2*.2);
680 G4double lp = std::log(p) - lmi;
681 G4double lp2 = lp*lp;
682 G4double rp2 = 1./p2;
683 elasticXsc = LE + ( pbe*lp2 + 6.72+32.6/p)/( 1. + rp2/p);
684 totalXsc = LE + ( pbt*lp2 + 38.2+52.7*rp2)/( 1. + 2.72*rp2*rp2);
685 }
686 }
687 else if( I==1 ) // np/pn
688 {
689 if( p < pmi )
690 {
691 G4double p2 = p*p;
692 elasticXsc = 1./( .00012 + p2*( .051 + .1*p2));
693 totalXsc = elasticXsc;
694 }
695 else if( p > pma )
696 {
697 G4double lp = std::log(p) - lmi;
698 G4double lp2 = lp*lp;
699 elasticXsc = pbe*lp2 + 6.72;
700 totalXsc = pbt*lp2 + 38.2;
701 }
702 else
703 {
704 G4double p2 = p*p;
705 G4double LE = 1./( .00012 + p2*( .051 + .1*p2 ) );
706 G4double lp = std::log(p) - lmi;
707 G4double lp2 = lp*lp;
708 G4double rp2 = 1./p2;
709 elasticXsc = LE + (pbe*lp2 + 6.72 + 30./p)/( 1. + .49*rp2/p);
710 totalXsc = LE + (pbt*lp2 + 38.2)/( 1. + .54*rp2*rp2);
711 }
712 }
713 else
714 {
715 G4cout<<"PDG incoding = "<<I<<" is not defined (0-1)"<<G4endl;
716
717 }
718 if( elasticXsc > totalXsc ) elasticXsc = totalXsc;
719
720 totalXsc *= millibarn;
721 elasticXsc *= millibarn;
722 inelasticXsc = totalXsc - elasticXsc;
723 if (inelasticXsc < 0.) inelasticXsc = 0.;
724
725 return inelasticXsc;
726}
727
728////////////////////////////////////////////////////////////////////////////////////
729//
730//
731
732G4double
733G4GGNuclNuclCrossSection::GetNucleusRadius( const G4DynamicParticle* ,
734 const G4Element* anElement)
735{
736 G4double At = anElement->GetN();
737 G4double oneThird = 1.0/3.0;
738 G4double cubicrAt = std::pow (At, oneThird);
739
740
741 G4double R; // = fRadiusConst*cubicrAt;
742 /*
743 G4double tmp = std::pow( cubicrAt-1., 3.);
744 tmp += At;
745 tmp *= 0.5;
746
747 if (At > 20.) // 20.
748 {
749 R = fRadiusConst*std::pow (tmp, oneThird);
750 }
751 else
752 {
753 R = fRadiusConst*cubicrAt;
754 }
755 */
756
757 R = fRadiusConst*cubicrAt;
758
759 // return R; // !!!!
760
761
762
763 G4double meanA = 21.;
764
765 G4double tauA1 = 40.;
766 G4double tauA2 = 10.;
767 G4double tauA3 = 5.;
768
769 G4double a1 = 0.85;
770 G4double b1 = 1. - a1;
771
772 G4double b2 = 0.3;
773 G4double b3 = 4.;
774
775 if (At > 20.) // 20.
776 {
777 R *= ( a1 + b1*std::exp( -(At - meanA)/tauA1) );
778 }
779 else if (At > 3.5)
780 {
781 R *= ( 1.0 + b2*( 1. - std::exp( (At - meanA)/tauA2) ) );
782 }
783 else
784 {
785 R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
786 }
787 return R;
788
789}
790
791////////////////////////////////////////////////////////////////////////////////////
792//
793//
794
795G4double
796G4GGNuclNuclCrossSection::GetNucleusRadius(G4double At)
797{
798 G4double R;
799
800 // R = GetNucleusRadiusGG(At);
801
802 R = GetNucleusRadiusDE(At);
803
804 return R;
805}
806
807///////////////////////////////////////////////////////////////////
808
809G4double
810G4GGNuclNuclCrossSection::GetNucleusRadiusGG(G4double At)
811{
812
813 G4double oneThird = 1.0/3.0;
814 G4double cubicrAt = std::pow (At, oneThird);
815
816
817 G4double R; // = fRadiusConst*cubicrAt;
818
819 /*
820 G4double tmp = std::pow( cubicrAt-1., 3.);
821 tmp += At;
822 tmp *= 0.5;
823
824 if (At > 20.)
825 {
826 R = fRadiusConst*std::pow (tmp, oneThird);
827 }
828 else
829 {
830 R = fRadiusConst*cubicrAt;
831 }
832 */
833
834 R = fRadiusConst*cubicrAt;
835
836 G4double meanA = 20.;
837 G4double tauA = 20.;
838
839 if ( At > 20.) // 20.
840 {
841 R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
842 }
843 else
844 {
845 R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
846 }
847
848 return R;
849
850}
851
852
853G4double
854G4GGNuclNuclCrossSection::GetNucleusRadiusDE(G4double A)
855{
856
857 // algorithm from diffuse-elastic
858
859 G4double R, r0, a11, a12, a13, a2, a3;
860
861 a11 = 1.26; // 1.08, 1.16
862 a12 = 1.; // 1.08, 1.16
863 a13 = 1.12; // 1.08, 1.16
864 a2 = 1.1;
865 a3 = 1.;
866
867
868 if( A < 50. )
869 {
870 if( 10 < A && A <= 15. ) r0 = a11*( 1 - std::pow(A, -2./3.) )*fermi; // 1.08*fermi;
871 else if( 15 < A && A <= 20 ) r0 = a12*( 1 - std::pow(A, -2./3.) )*fermi;
872 else if( 20 < A && A <= 30 ) r0 = a13*( 1 - std::pow(A, -2./3.) )*fermi;
873 else r0 = a2*fermi;
874
875 R = r0*std::pow( A, 1./3. );
876 }
877 else
878 {
879 r0 = a3*fermi;
880
881 R = r0*std::pow(A, 0.27);
882 }
883 return R;
884
885
886
887}
888
889
890
891
892
893
894////////////////////////////////////////////////////////////////////////////////////
895//
896//
897
898G4double G4GGNuclNuclCrossSection::CalculateEcmValue( const G4double mp ,
899 const G4double mt ,
900 const G4double Plab )
901{
902 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
903 G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
904 // G4double Pcm = Plab * mt / Ecm;
905 // G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
906
907 return Ecm ; // KEcm;
908}
909
910
911////////////////////////////////////////////////////////////////////////////////////
912//
913//
914
915G4double G4GGNuclNuclCrossSection::CalcMandelstamS( const G4double mp ,
916 const G4double mt ,
917 const G4double Plab )
918{
919 G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
920 G4double sMand = mp*mp + mt*mt + 2*Elab*mt ;
921
922 return sMand;
923}
924
925
926//
927//
928///////////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.