58 FTFhNcmsEnergy( 0.0 ),
60 FTFXtotal( 0.0 ), FTFXelastic( 0.0 ), FTFXinelastic( 0.0 ), FTFXannihilation( 0.0 ),
61 ProbabilityOfAnnihilation( 0.0 ), ProbabilityOfElasticScatt( 0.0 ),
62 RadiusOfHNinteractions2( 0.0 ), FTFSlope( 0.0 ),
63 AvaragePt2ofElasticScattering( 0.0 ), FTFGamma0( 0.0 ),
64 DeltaProbAtQuarkExchange( 0.0 ), ProbOfSameQuarkExchange( 0.0 ),
65 ProjMinDiffMass( 0.0 ), ProjMinNonDiffMass( 0.0 ),
66 TarMinDiffMass( 0.0 ), TarMinNonDiffMass( 0.0 ),
67 AveragePt2( 0.0 ), ProbLogDistr( 0.0 ),
69 MaxNumberOfCollisions( 0.0 ), ProbOfInelInteraction( 0.0 ), CofNuclearDestruction( 0.0 ),
70 R2ofNuclearDestruction( 0.0 ), ExcitationEnergyPerWoundedNucleon( 0.0 ),
71 DofNuclearDestruction( 0.0 ), Pt2ofNuclearDestruction( 0.0 ), MaxPt2ofNuclearDestruction( 0.0 )
73 for (
G4int i = 0; i < 4; i++ ) {
74 for (
G4int j = 0; j < 7; j++ ) {
88 G4ThreadLocal bool G4FTFParameters::chipsComponentXSisInitialized =
false;
102 G4int ProjectileabsPDGcode = std::abs( ProjectilePDGcode );
104 G4double ProjectileMass2 = ProjectileMass * ProjectileMass;
106 G4int ProjectileBaryonNumber( 0 ), AbsProjectileBaryonNumber( 0 ), AbsProjectileCharge( 0 );
107 G4bool ProjectileIsNucleus =
false;
110 ProjectileIsNucleus =
true;
112 AbsProjectileBaryonNumber = std::abs( ProjectileBaryonNumber );
114 if ( ProjectileBaryonNumber > 1 ) {
115 ProjectilePDGcode = 2212; ProjectileabsPDGcode = 2212;
117 ProjectilePDGcode = -2212; ProjectileabsPDGcode = 2212;
120 ProjectileMass2 =
sqr( ProjectileMass );
124 G4double TargetMass2 = TargetMass * TargetMass;
127 G4double Elab = std::sqrt( Plab*Plab + ProjectileMass2 );
128 G4double KineticEnergy = Elab - ProjectileMass;
130 G4double S = ProjectileMass2 + TargetMass2 + 2.0*TargetMass*Elab;
132 #ifdef debugFTFparams
133 G4cout <<
"--------- FTF Parameters --------------" <<
G4endl <<
"Proj Plab "
134 << ProjectilePDGcode <<
" " << Plab <<
G4endl <<
"Mass KinE " << ProjectileMass
135 <<
" " << KineticEnergy <<
G4endl <<
" A Z " << theA <<
" " << theZ <<
G4endl;
138 G4double Ylab, Xtotal, Xelastic, Xannihilation;
139 G4int NumberOfTargetNucleons;
141 Ylab = 0.5 * std::log( (Elab + Plab)/(Elab - Plab) );
146 #ifdef debugFTFparams
150 TargetMass /=
GeV; TargetMass2 /= (
GeV*
GeV);
151 ProjectileMass /=
GeV; ProjectileMass2 /= (
GeV*
GeV);
157 if ( ! chipsComponentXSisInitialized ) {
158 chipsComponentXSisInitialized =
true;
169 G4int NumberOfTargetProtons = theZ;
170 G4int NumberOfTargetNeutrons = theA - theZ;
171 NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
173 if ( ProjectilePDGcode == 2212 || ProjectilePDGcode == 2112 ) {
181 #ifdef debugFTFparams
186 if ( ! ProjectileIsNucleus ) {
187 Xtotal = ( NumberOfTargetProtons * XtotPP + NumberOfTargetNeutrons * XtotPN ) /
188 NumberOfTargetNucleons;
189 Xelastic = ( NumberOfTargetProtons * XelPP + NumberOfTargetNeutrons * XelPN ) /
190 NumberOfTargetNucleons;
193 AbsProjectileCharge * NumberOfTargetProtons * XtotPP +
194 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
195 NumberOfTargetNeutrons * XtotPP
197 ( AbsProjectileCharge * NumberOfTargetNeutrons +
198 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
199 NumberOfTargetProtons ) * XtotPN
200 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
202 AbsProjectileCharge * NumberOfTargetProtons * XelPP +
203 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
204 NumberOfTargetNeutrons * XelPP
206 ( AbsProjectileCharge * NumberOfTargetNeutrons +
207 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
208 NumberOfTargetProtons ) * XelPN
209 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
216 }
else if ( ProjectilePDGcode < -1000 ) {
218 G4double X_a( 0.0 ), X_b( 0.0 ), X_c( 0.0 ), X_d( 0.0 );
219 G4double MesonProdThreshold = ProjectileMass + TargetMass +
220 ( 2.0 * 0.14 + 0.016 );
222 if ( PlabPerParticle < 40.0*
MeV ) {
230 G4double LogS = std::log( ECMSsqr / 33.0625 );
231 G4double Xasmpt = 36.04 + 0.304*LogS*LogS;
232 LogS = std::log( SqrtS / 20.74 );
233 G4double Basmpt = 11.92 + 0.3036*LogS*LogS;
234 G4double R0 = std::sqrt( 0.40874044*Xasmpt - Basmpt );
236 G4double FlowF = SqrtS / std::sqrt( ECMSsqr*ECMSsqr + ProjectileMass2*ProjectileMass2 +
237 TargetMass2*TargetMass2 - 2.0*ECMSsqr*ProjectileMass2
238 - 2.0*ECMSsqr*TargetMass2
239 - 2.0*ProjectileMass2*TargetMass2 );
241 Xtotal = Xasmpt * ( 1.0 + 13.55*FlowF/R0/R0/R0*
242 (1.0 - 4.47/SqrtS + 12.38/ECMSsqr - 12.43/SqrtS/ECMSsqr) );
244 Xasmpt = 4.4 + 0.101*LogS*LogS;
245 Xelastic = Xasmpt * ( 1.0 + 59.27*FlowF/R0/R0/R0*
246 (1.0 - 6.95/SqrtS + 23.54/ECMSsqr - 25.34/SqrtS/ECMSsqr ) );
255 if ( SqrtS < MesonProdThreshold ) {
256 X_b = 3.13 + 140.0*std::pow( MesonProdThreshold - SqrtS, 2.5 );
263 X_c = 2.0*FlowF*
sqr( ProjectileMass + TargetMass )/ECMSsqr;
275 G4double Xann_on_P( 0.0), Xann_on_N( 0.0 );
277 if ( ProjectilePDGcode == -2212 ) {
278 Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0;
279 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
280 }
else if ( ProjectilePDGcode == -2112 ) {
281 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
282 Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0;
283 }
else if ( ProjectilePDGcode == -3122 ) {
284 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
285 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
286 }
else if ( ProjectilePDGcode == -3112 ) {
287 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
288 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
289 }
else if ( ProjectilePDGcode == -3212 ) {
290 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
291 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
292 }
else if ( ProjectilePDGcode == -3222 ) {
293 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
294 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
295 }
else if ( ProjectilePDGcode == -3312 ) {
296 Xann_on_P = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
297 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
298 }
else if ( ProjectilePDGcode == -3322 ) {
299 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
300 Xann_on_N = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
301 }
else if ( ProjectilePDGcode == -3334 ) {
302 Xann_on_P = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
303 Xann_on_N = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
305 G4cout <<
"Unknown anti-baryon for FTF annihilation" <<
G4endl;
310 if ( ! ProjectileIsNucleus ) {
311 Xannihilation = ( NumberOfTargetProtons * Xann_on_P + NumberOfTargetNeutrons * Xann_on_N )
312 / NumberOfTargetNucleons;
315 ( AbsProjectileCharge * NumberOfTargetProtons +
316 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
317 NumberOfTargetNeutrons ) * Xann_on_P
319 ( AbsProjectileCharge * NumberOfTargetNeutrons +
320 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
321 NumberOfTargetProtons ) * Xann_on_N
322 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
326 MesonProdThreshold = ProjectileMass + TargetMass + (0.14 + 0.08);
327 if ( SqrtS > MesonProdThreshold ) {
328 Xftf = 36.0 * ( 1.0 - MesonProdThreshold/SqrtS );
331 Xtotal = Xelastic + Xannihilation + Xftf;
333 #ifdef debugFTFparams
334 G4cout <<
"Plab Xtotal, Xelastic Xinel Xftf " << Plab <<
" " << Xtotal <<
" " << Xelastic
335 <<
" " << Xtotal - Xelastic <<
" " << Xtotal - Xelastic - Xannihilation << G4endl
336 <<
"Plab Xelastic/Xtotal, Xann/Xin " << Plab <<
" " << Xelastic/Xtotal <<
" "
337 << Xannihilation/(Xtotal - Xelastic) << G4endl;
340 }
else if ( ProjectilePDGcode == 211 ) {
347 Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
348 / NumberOfTargetNucleons;
349 Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
350 / NumberOfTargetNucleons;
355 }
else if ( ProjectilePDGcode == -211 ) {
362 Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
363 / NumberOfTargetNucleons;
364 Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
365 / NumberOfTargetNucleons;
370 }
else if ( ProjectilePDGcode == 111 ) {
378 G4double XtotPiP = ( XtotPipP + XtotPimP ) / 2.0;
380 G4double XelPiP = ( XelPipP + XelPimP ) / 2.0;
382 Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
383 / NumberOfTargetNucleons;
384 Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
385 / NumberOfTargetNucleons;
390 }
else if ( ProjectilePDGcode == 321 ) {
397 Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
398 / NumberOfTargetNucleons;
399 Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
400 / NumberOfTargetNucleons;
405 }
else if ( ProjectilePDGcode == -321 ) {
412 Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
413 / NumberOfTargetNucleons;
414 Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
415 / NumberOfTargetNucleons;
420 }
else if ( ProjectilePDGcode == 311 || ProjectilePDGcode == 130 ||
421 ProjectilePDGcode == 310 ) {
429 G4double XtotKP = ( XtotKpP + XtotKmP ) / 2.0;
431 G4double XelKP = ( XelKpP + XelKmP ) / 2.0;
433 Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
434 / NumberOfTargetNucleons;
435 Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
436 / NumberOfTargetNucleons;
449 Xtotal = ( NumberOfTargetProtons * XtotPP + NumberOfTargetNeutrons * XtotPN )
450 / NumberOfTargetNucleons;
451 Xelastic = ( NumberOfTargetProtons * XelPP + NumberOfTargetNeutrons * XelPN )
452 / NumberOfTargetNucleons;
478 if ( Xtotal - Xelastic == 0.0 ) {
490 SetSlope( Xtotal*Xtotal/16.0/
pi/Xelastic/0.3894 );
505 if ( ProjectilePDGcode > 1000 ) {
508 SetParams( 0, 13.71, 1.75, -214.4 , 4.25, 0.0, 0.632, 1.45 );
509 SetParams( 1, 0.2, 0.0 , - 3.289, 2.0 , 0.0, 0.0 , 1.40 );
510 SetParams( 2, 6.0/Xinel, 0.0 , -6.0/Xinel*8.48, 2.25, 0.0, 0.0 , 0.95 );
511 SetParams( 3, 6.0/Xinel, 0.0 , -6.0/Xinel*8.48, 2.25, 0.0, 0.0 , 0.95 );
513 SetParams( 1, 0.3, 0.5 , - 3.289, 2.0 , 0.0, 0.0 , 1.40 );
514 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
515 SetParams( 2, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
516 SetParams( 3, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
519 if ( NumberOfTargetNucleons > 26 ) {
531 }
else if( ProjectilePDGcode < -1000 ) {
534 SetParams( 0, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 );
535 SetParams( 1, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 );
537 SetParams( 2, 6.0/Xftf, 0.0 , -6.0/Xftf*8.48, 2.25, 0.0, 0.0 , 0.95 );
538 SetParams( 3, 6.0/Xftf, 0.0 , -6.0/Xftf*8.48, 2.25, 0.0, 0.0 , 0.95 );
540 SetParams( 2, 0.5 , 0.0 , 0.0 , 0.0 , 0.0, 0.5 , 1000.0 );
541 SetParams( 3, 0.5 , 0.0 , 0.0 , 0.0 , 0.0, 0.5 , 1000.0 );
543 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
544 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
545 SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
556 }
else if ( ProjectileabsPDGcode == 211 || ProjectilePDGcode == 111 ) {
559 SetParams( 0, 568.0 , 2.1 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
560 SetParams( 1, 6.0 , 0.6 , -26.9 , 1.1 , 0.0, 0.0 , 3.0 );
562 if ( Xinel > 0.0 ) Wprd = 0.64 *( 6.2 - 3.7*std::exp( -
sqr( SqrtS - 7.0 ) / 16.0 ) ) / Xinel;
563 SetParams( 2, Wprd, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
565 if ( Xinel > 0.0 ) Wtrd = ( 2.0 + 22.0/ECMSsqr ) / Xinel;
566 SetParams( 3, Wtrd, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
567 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
568 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
569 SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
579 }
else if ( ProjectileabsPDGcode == 321 || ProjectileabsPDGcode == 311 ||
580 ProjectilePDGcode == 130 || ProjectilePDGcode == 310 ) {
583 SetParams( 0, 70.0 , 2.75, 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
584 SetParams( 1, 30.0 , 1.5 , -50.57 , 1.83, 0.0, 0.0 , 1.70 );
585 SetParams( 2, 0.6 , 0.75, -12.05 , 3.25, 0.0, 0.0 , 1.20 );
586 SetParams( 3, 6.0 , 1.0 , -12.08 , 1.5 , 0.1, 0.0 , 1.20 );
587 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
588 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
589 SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
602 SetParams( 0, 13.71, 1.75, -214.4 , 4.25, 0.0, 0.632, 1.45 );
603 SetParams( 1, 0.2 , 0.0 , -16.445, 2.0 , 0.0, 0.0 , 1.40 );
604 SetParams( 2, 6.0/Xinel, 0.0 , -6.0/Xinel*8.48, 2.25, 0.0, 0.0 , 0.95 );
605 SetParams( 3, 6.0/Xinel, 0.0 , -6.0/Xinel*8.48, 2.25, 0.0, 0.0 , 0.95 );
606 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
607 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
608 SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 );
626 G4double Puubar( 1.0/3.0 ), Pddbar( 1.0/3.0 ), Pssbar( 1.0/3.0 );
631 if ( ProjectileabsPDGcode < 1000 ) {
634 ( 1.0 + std::exp( 4.0*(Ylab - 2.1) ) ) );
638 ( 1.0 + std::exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
641 }
else if ( ProjectilePDGcode < -1000 ) {
645 ( 1.0 + std::exp( 4.0*(Ylab - 2.1) ) ) );
649 ( 1.0 + std::exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
664 ( 1.0 + std::exp( 4.0*(Ylab - 2.1) ) ) );
668 ( 1.0 + std::exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
void SetProbLogDistr(const G4double aValue)
G4int GetPDGEncoding() const
void SetTarMinNonDiffMass(const G4double aValue)
void SetDofNuclearDestruction(const G4double aValue)
void SetSlope(const G4double Slope)
G4ChipsComponentXS * FTFxsManager
static G4KaonMinus * KaonMinus()
void SetGamma0(const G4double Gamma0)
G4GLOB_DLL std::ostream G4cout
void SetAveragePt2(const G4double aValue)
G4double ProcParams[4][7]
void SetElastisCrossSection(const G4double Xelastic)
void SetPt2Kink(const G4double aValue)
void SetRadiusOfHNinteractions2(const G4double Radius2)
void SetMaxNumberOfCollisions(const G4double aValue, const G4double bValue)
void SetExcitationEnergyPerWoundedNucleon(const G4double aValue)
static G4Proton * Proton()
void SetProjMinNonDiffMass(const G4double aValue)
static G4PionPlus * PionPlus()
void SetQuarkProbabilitiesAtGluonSplitUp(const G4double Puubar, const G4double Pddbar, const G4double Pssbar)
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double N)
static G4Neutron * Neutron()
G4double GetProcProb(const G4int ProcN, const G4double y)
void SetInelasticCrossSection(const G4double Xinelastic)
void SetProbabilityOfElasticScatt(const G4double Xtotal, const G4double Xelastic)
void SetDeltaProbAtQuarkExchange(const G4double aValue)
void SetTotalCrossSection(const G4double Xtotal)
void SetProjMinDiffMass(const G4double aValue)
G4double ProbOfSameQuarkExchange
G4double GetPDGMass() const
void SetMaxPt2ofNuclearDestruction(const G4double aValue)
void SetParams(const G4int ProcN, const G4double A1, const G4double B1, const G4double A2, const G4double B2, const G4double A3, const G4double Atop, const G4double Ymin)
static G4PionMinus * PionMinus()
void SetCofNuclearDestruction(const G4double aValue)
void SetProbOfSameQuarkExchange(const G4double aValue)
G4double FTFXannihilation
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double N)
void SetProbabilityOfAnnihilation(const G4double aValue)
static G4KaonPlus * KaonPlus()
G4double GetPDGCharge() const
void SetPt2ofNuclearDestruction(const G4double aValue)
void SetAvaragePt2ofElasticScattering(const G4double aPt2)
void SetTarMinDiffMass(const G4double aValue)
void SetR2ofNuclearDestruction(const G4double aValue)
G4int GetBaryonNumber() const