Geant4  10.01.p01
G4FTFParameters.cc
Go to the documentation of this file.
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 // $Id: G4FTFParameters.cc 86868 2014-11-19 14:46:25Z gcosmo $
28 // GEANT4 tag $Name: $
29 //
30 
31 #include <utility>
32 
33 #include "G4FTFParameters.hh"
34 
35 #include "G4ios.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 
39 #include "G4ParticleDefinition.hh"
40 
41 #include "G4Proton.hh"
42 #include "G4Neutron.hh"
43 
44 #include "G4PionPlus.hh"
45 #include "G4PionMinus.hh"
46 #include "G4KaonPlus.hh"
47 #include "G4KaonMinus.hh"
48 
49 
50 //============================================================================
51 
52 //#define debugFTFparams
53 
54 
55 //============================================================================
56 
58  FTFhNcmsEnergy( 0.0 ),
59  FTFxsManager( 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 ), ProbLogDistrPrD(0.0), // Uzhi Oct 2014
66  TarMinDiffMass( 0.0 ), TarMinNonDiffMass( 0.0 ),
67  AveragePt2( 0.0 ), ProbLogDistr( 0.0 ),
68  Pt2kink( 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 )
72 {
73  for ( G4int i = 0; i < 4; i++ ) {
74  for ( G4int j = 0; j < 7; j++ ) {
75  ProcParams[i][j] = 0.0;
76  }
77  }
78 }
79 
80 
81 //============================================================================
82 
84 
85 
86 //============================================================================
87 
90 
91 
92 //============================================================================
93 
95  G4int theA, G4int theZ, G4double PlabPerParticle ) {
96 
97  FTFXannihilation = 0.0;
98  FTFhNcmsEnergy = 0.0;
100 
101  G4int ProjectilePDGcode = particle->GetPDGEncoding();
102  G4int ProjectileabsPDGcode = std::abs( ProjectilePDGcode );
103  G4double ProjectileMass = particle->GetPDGMass();
104  G4double ProjectileMass2 = ProjectileMass * ProjectileMass;
105 
106  G4int ProjectileBaryonNumber( 0 ), AbsProjectileBaryonNumber( 0 ), AbsProjectileCharge( 0 );
107  G4bool ProjectileIsNucleus = false;
108 
109  if ( std::abs( particle->GetBaryonNumber() ) > 1 ) { // The projectile is a nucleus
110  ProjectileIsNucleus = true;
111  ProjectileBaryonNumber = particle->GetBaryonNumber();
112  AbsProjectileBaryonNumber = std::abs( ProjectileBaryonNumber );
113  AbsProjectileCharge = G4int( particle->GetPDGCharge() );
114  if ( ProjectileBaryonNumber > 1 ) {
115  ProjectilePDGcode = 2212; ProjectileabsPDGcode = 2212; // Proton
116  } else {
117  ProjectilePDGcode = -2212; ProjectileabsPDGcode = 2212; // Anti-Proton
118  }
119  ProjectileMass = G4Proton::Proton()->GetPDGMass();
120  ProjectileMass2 = sqr( ProjectileMass );
121  }
122 
123  G4double TargetMass = G4Proton::Proton()->GetPDGMass();
124  G4double TargetMass2 = TargetMass * TargetMass;
125 
126  G4double Plab = PlabPerParticle;
127  G4double Elab = std::sqrt( Plab*Plab + ProjectileMass2 );
128  G4double KineticEnergy = Elab - ProjectileMass;
129 
130  G4double S = ProjectileMass2 + TargetMass2 + 2.0*TargetMass*Elab;
131 
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;
136  #endif
137 
138  G4double Ylab, Xtotal, Xelastic, Xannihilation;
139  G4int NumberOfTargetNucleons;
140 
141  Ylab = 0.5 * std::log( (Elab + Plab)/(Elab - Plab) );
142 
143  G4double ECMSsqr = S/GeV/GeV;
144  G4double SqrtS = std::sqrt( S )/GeV;
145 
146  #ifdef debugFTFparams
147  G4cout << "Sqrt(s) " << SqrtS << G4endl;
148  #endif
149 
150  TargetMass /= GeV; TargetMass2 /= (GeV*GeV);
151  ProjectileMass /= GeV; ProjectileMass2 /= (GeV*GeV);
152 
153  // Andrea Dotti (13Jan2013):
154  // The following lines are changed for G4MT. Originally the code was:
155  // static G4ChipsComponentXS* _instance = new G4ChipsComponentXS(); // Witek Pokorski
156  // Note the code could go back at original if _instance could be shared among threads
160  }
162  FTFxsManager = _instance;
163 
164  Plab /= GeV;
165  G4double Xftf = 0.0;
166  //G4double LogPlab = std::log( Plab );
167  //G4double sqrLogPlab = LogPlab * LogPlab;
168 
169  G4int NumberOfTargetProtons = theZ;
170  G4int NumberOfTargetNeutrons = theA - theZ;
171  NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
172 
173  if ( ProjectilePDGcode == 2212 || ProjectilePDGcode == 2112 ) { // Projectile is nucleon
174 
175  G4double XtotPP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
177  G4double XtotPN = FTFxsManager->GetTotalElementCrossSection( Neutron, KineticEnergy, 1, 0 );
178  G4double XelPP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
179  G4double XelPN = FTFxsManager->GetElasticElementCrossSection( Neutron, KineticEnergy, 1, 0 );
180 
181  #ifdef debugFTFparams
182  G4cout << "XsPP " << XtotPP/millibarn << " " << XelPP/millibarn << G4endl
183  << "XsPN " << XtotPN/millibarn << " " << XelPN/millibarn << G4endl;
184  #endif
185 
186  if ( ! ProjectileIsNucleus ) { // Projectile is hadron
187  Xtotal = ( NumberOfTargetProtons * XtotPP + NumberOfTargetNeutrons * XtotPN ) /
188  NumberOfTargetNucleons;
189  Xelastic = ( NumberOfTargetProtons * XelPP + NumberOfTargetNeutrons * XelPN ) /
190  NumberOfTargetNucleons;
191  } else { // Projectile is a nucleus
192  Xtotal = (
193  AbsProjectileCharge * NumberOfTargetProtons * XtotPP +
194  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
195  NumberOfTargetNeutrons * XtotPP
196  +
197  ( AbsProjectileCharge * NumberOfTargetNeutrons +
198  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
199  NumberOfTargetProtons ) * XtotPN
200  ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
201  Xelastic= (
202  AbsProjectileCharge * NumberOfTargetProtons * XelPP +
203  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
204  NumberOfTargetNeutrons * XelPP
205  +
206  ( AbsProjectileCharge * NumberOfTargetNeutrons +
207  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
208  NumberOfTargetProtons ) * XelPN
209  ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
210  }
211 
212  Xannihilation = 0.0;
213  Xtotal /= millibarn;
214  Xelastic /= millibarn;
215 
216  } else if ( ProjectilePDGcode < -1000 ) { // Projectile is anti_baryon
217 
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 ); // 2 Mpi + DeltaE;
221 
222  if ( PlabPerParticle < 40.0*MeV ) { // Low energy limits. Projectile at rest.
223  Xtotal = 1512.9; // mb
224  Xelastic = 473.2; // mb
225  X_a = 625.1; // mb
226  X_b = 9.780; // mb
227  X_c = 49.989; // mb
228  X_d = 6.614; // mb
229  } else { // Total and elastic cross section of PbarP interactions a'la Arkhipov
230  G4double LogS = std::log( ECMSsqr / 33.0625 );
231  G4double Xasmpt = 36.04 + 0.304*LogS*LogS; // mb
232  LogS = std::log( SqrtS / 20.74 );
233  G4double Basmpt = 11.92 + 0.3036*LogS*LogS; // GeV^(-2)
234  G4double R0 = std::sqrt( 0.40874044*Xasmpt - Basmpt ); // GeV^(-1)
235 
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 );
240 
241  Xtotal = Xasmpt * ( 1.0 + 13.55*FlowF/R0/R0/R0*
242  (1.0 - 4.47/SqrtS + 12.38/ECMSsqr - 12.43/SqrtS/ECMSsqr) ); // mb
243 
244  Xasmpt = 4.4 + 0.101*LogS*LogS; // mb
245  Xelastic = Xasmpt * ( 1.0 + 59.27*FlowF/R0/R0/R0*
246  (1.0 - 6.95/SqrtS + 23.54/ECMSsqr - 25.34/SqrtS/ECMSsqr ) ); // mb
247 
248  //G4cout << "Param Xtotal Xelastic " << Xtotal << " " << Xelastic << G4endl
249  // << "FlowF " << FlowF << " SqrtS " << SqrtS << G4endl
250  // << "Param Xelastic-NaN " << Xelastic << " "
251  // << 1.5*16.654/pow(ECMSsqr/2.176/2.176,2.2) << " " << ECMSsqr << G4endl;
252 
253  X_a = 25.0*FlowF; // mb, 3-shirts diagram
254 
255  if ( SqrtS < MesonProdThreshold ) {
256  X_b = 3.13 + 140.0*std::pow( MesonProdThreshold - SqrtS, 2.5 ); // mb anti-quark-quark annihilation
257  Xelastic -= 3.0*X_b; // Xel-X(PbarP->NNbar)
258  } else {
259  X_b = 6.8/SqrtS; // mb anti-quark-quark annihilation
260  Xelastic -= 3.0*X_b; // Xel-X(PbarP->NNbar)
261  }
262 
263  X_c = 2.0*FlowF*sqr( ProjectileMass + TargetMass )/ECMSsqr; // mb rearrangement
264 
265  //G4cout << "Old new Xa " << 35.*FlowF << " " << 25.*FlowF << G4endl;
266 
267  X_d = 23.3/ECMSsqr; // mb anti-quark-quark string creation
268  }
269 
270  //G4cout << "Param Xtotal Xelastic " << Xtotal << " " << Xelastic << G4endl
271  // << "Para a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d << G4endl;
272  // << "Para a b c d " << X_a << " " << 5.*X_b << " " << 5.*X_c << " " << 6.*X_d
273  // << G4endl;
274 
275  G4double Xann_on_P( 0.0), Xann_on_N( 0.0 );
276 
277  if ( ProjectilePDGcode == -2212 ) { // Pbar+P/N
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 ) { // NeutrBar+P/N
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 ) { // LambdaBar+P/N
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 ) { // Sigma-Bar+P/N
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 ) { // Sigma0Bar+P/N
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 ) { // Sigma+Bar+P/N
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 ) { // Xi-Bar+P/N
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 ) { // Xi0Bar+P/N
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 ) { // Omega-Bar+P/N
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;
304  } else {
305  G4cout << "Unknown anti-baryon for FTF annihilation" << G4endl;
306  }
307 
308  //G4cout << "Sum " << Xann_on_P << G4endl;
309 
310  if ( ! ProjectileIsNucleus ) { // Projectile is anti-baryon
311  Xannihilation = ( NumberOfTargetProtons * Xann_on_P + NumberOfTargetNeutrons * Xann_on_N )
312  / NumberOfTargetNucleons;
313  } else { // Projectile is a nucleus
314  Xannihilation = (
315  ( AbsProjectileCharge * NumberOfTargetProtons +
316  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
317  NumberOfTargetNeutrons ) * Xann_on_P
318  +
319  ( AbsProjectileCharge * NumberOfTargetNeutrons +
320  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
321  NumberOfTargetProtons ) * Xann_on_N
322  ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
323  }
324 
325  //G4double Xftf = 0.0;
326  MesonProdThreshold = ProjectileMass + TargetMass + (0.14 + 0.08); // Mpi + DeltaE
327  if ( SqrtS > MesonProdThreshold ) {
328  Xftf = 36.0 * ( 1.0 - MesonProdThreshold/SqrtS );
329  }
330 
331  Xtotal = Xelastic + Xannihilation + Xftf;
332 
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;
338  #endif
339 
340  } else if ( ProjectilePDGcode == 211 ) { // Projectile is PionPlus
341 
342  G4double XtotPiP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
344  G4double XtotPiN = FTFxsManager->GetTotalElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
345  G4double XelPiP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
346  G4double XelPiN = FTFxsManager->GetElasticElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
347  Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
348  / NumberOfTargetNucleons;
349  Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
350  / NumberOfTargetNucleons;
351  Xannihilation = 0.0;
352  Xtotal /= millibarn;
353  Xelastic /= millibarn;
354 
355  } else if ( ProjectilePDGcode == -211 ) { // Projectile is PionMinus
356 
357  G4double XtotPiP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
359  G4double XtotPiN = FTFxsManager->GetTotalElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
360  G4double XelPiP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
361  G4double XelPiN = FTFxsManager->GetElasticElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
362  Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
363  / NumberOfTargetNucleons;
364  Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
365  / NumberOfTargetNucleons;
366  Xannihilation = 0.0;
367  Xtotal /= millibarn;
368  Xelastic /= millibarn;
369 
370  } else if ( ProjectilePDGcode == 111 ) { // Projectile is PionZero
371 
373  G4double XtotPipP = FTFxsManager->GetTotalElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
375  G4double XtotPimP = FTFxsManager->GetTotalElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
376  G4double XelPipP = FTFxsManager->GetElasticElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
377  G4double XelPimP = FTFxsManager->GetElasticElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
378  G4double XtotPiP = ( XtotPipP + XtotPimP ) / 2.0;
379  G4double XtotPiN = XtotPiP;
380  G4double XelPiP = ( XelPipP + XelPimP ) / 2.0;
381  G4double XelPiN = XelPiP;
382  Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
383  / NumberOfTargetNucleons;
384  Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
385  / NumberOfTargetNucleons;
386  Xannihilation = 0.0;
387  Xtotal /= millibarn;
388  Xelastic /= millibarn;
389 
390  } else if ( ProjectilePDGcode == 321 ) { // Projectile is KaonPlus
391 
392  G4double XtotKP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
394  G4double XtotKN = FTFxsManager->GetTotalElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
395  G4double XelKP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
396  G4double XelKN = FTFxsManager->GetElasticElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
397  Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
398  / NumberOfTargetNucleons;
399  Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
400  / NumberOfTargetNucleons;
401  Xannihilation = 0.0;
402  Xtotal /= millibarn;
403  Xelastic /= millibarn;
404 
405  } else if ( ProjectilePDGcode == -321 ) { // Projectile is KaonMinus
406 
407  G4double XtotKP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
409  G4double XtotKN = FTFxsManager->GetTotalElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
410  G4double XelKP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
411  G4double XelKN = FTFxsManager->GetElasticElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
412  Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
413  / NumberOfTargetNucleons;
414  Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
415  / NumberOfTargetNucleons;
416  Xannihilation = 0.0;
417  Xtotal /= millibarn;
418  Xelastic /= millibarn;
419 
420  } else if ( ProjectilePDGcode == 311 || ProjectilePDGcode == 130 ||
421  ProjectilePDGcode == 310 ) { // Projectile is KaonZero
422 
424  G4double XtotKpP = FTFxsManager->GetTotalElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
426  G4double XtotKmP = FTFxsManager->GetTotalElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
427  G4double XelKpP = FTFxsManager->GetElasticElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
428  G4double XelKmP = FTFxsManager->GetElasticElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
429  G4double XtotKP = ( XtotKpP + XtotKmP ) / 2.0;
430  G4double XtotKN = XtotKP;
431  G4double XelKP = ( XelKpP + XelKmP ) / 2.0;
432  G4double XelKN = XelKP;
433  Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
434  / NumberOfTargetNucleons;
435  Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
436  / NumberOfTargetNucleons;
437  Xannihilation = 0.0;
438  Xtotal /= millibarn;
439  Xelastic /= millibarn;
440 
441  } else { // Projectile is undefined, Nucleon assumed
442 
444  G4double XtotPP = FTFxsManager->GetTotalElementCrossSection( Proton, KineticEnergy, 1, 0 );
446  G4double XtotPN = FTFxsManager->GetTotalElementCrossSection( Neutron, KineticEnergy, 1, 0 );
447  G4double XelPP = FTFxsManager->GetElasticElementCrossSection( Proton, KineticEnergy, 1, 0 );
448  G4double XelPN = FTFxsManager->GetElasticElementCrossSection( Neutron, KineticEnergy, 1, 0 );
449  Xtotal = ( NumberOfTargetProtons * XtotPP + NumberOfTargetNeutrons * XtotPN )
450  / NumberOfTargetNucleons;
451  Xelastic = ( NumberOfTargetProtons * XelPP + NumberOfTargetNeutrons * XelPN )
452  / NumberOfTargetNucleons;
453  Xannihilation = 0.0;
454  Xtotal /= millibarn;
455  Xelastic /= millibarn;
456 
457  };
458 
459  // Geometrical parameters
460  SetTotalCrossSection( Xtotal );
461  SetElastisCrossSection( Xelastic );
462  SetInelasticCrossSection( Xtotal - Xelastic );
463 
464  //G4cout << "Plab Xtotal, Xelastic Xinel Xftf " << Plab << " " << Xtotal << " " << Xelastic
465  // << " " << Xtotal - Xelastic << " " << Xtotal - Xelastic - Xannihilation << G4endl;
466  //if (Xtotal - Xelastic != 0.0 ) {
467  // G4cout << "Plab Xelastic/Xtotal, Xann/Xin " << Plab << " " << Xelastic/Xtotal
468  // << " " << Xannihilation / (Xtotal - Xelastic) << G4endl;
469  //} else {
470  // G4cout << "Plab Xelastic/Xtotal, Xann " << Plab << " " << Xelastic/Xtotal
471  // << " " << Xannihilation << G4endl;
472  //}
473  //G4int Uzhi; G4cin >> Uzhi;
474 
475  // Interactions with elastic and inelastic collisions
476  SetProbabilityOfElasticScatt( Xtotal, Xelastic );
477  SetRadiusOfHNinteractions2( Xtotal/pi/10.0 );
478  if ( Xtotal - Xelastic == 0.0 ) {
480  } else {
481  SetProbabilityOfAnnihilation( Xannihilation / (Xtotal - Xelastic) );
482  }
483 
484  // No elastic scattering
485  //SetProbabilityOfElasticScatt( Xtotal, 0.0 );
486  //SetRadiusOfHNinteractions2( (Xtotal - Xelastic)/pi/10.0 );
487  //SetProbabilityOfAnnihilation( 1.0 );
488  //SetProbabilityOfAnnihilation( 0.0 );
489 
490  SetSlope( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
491  // (GeV/c)^(-2))
492  //G4cout << "Slope " << GetSlope() << G4endl;
493  SetGamma0( GetSlope()*Xtotal/10.0/2.0/pi );
494 
495  // Parameters of elastic scattering
496  // Gaussian parametrization of elastic scattering amplitude assumed
497  SetAvaragePt2ofElasticScattering( 1.0/( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 )*GeV*GeV );
498  //G4cout << "AvaragePt2ofElasticScattering " << GetAvaragePt2ofElasticScattering() << G4endl;
499 
500  // Parameters of excitations
501 
502  G4double Xinel = Xtotal - Xelastic; // Uzhi 25.04.2012
503  //G4cout << "Param ProjectilePDGcode " << ProjectilePDGcode << G4endl;
504 
505  if ( ProjectilePDGcode > 1000 ) { // Projectile is baryon
506  // Proc# A1 B1 A2 B2 A3 Atop Ymin
507  SetParams( 0, 13.71, 1.75, -214.5, 4.25, 0.0, 0.5 , 1.1 ); // Qexchange without Exc.
508  SetParams( 1, 25.0, 1.0, -50.34, 1.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc.
509  SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Projectile diffraction
510  SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Target diffraction
511  SetParams( 4, 1.0, 0.0 , -2.01 , 0.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc. Additional multiply
512 //
513  if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
514  SetParams( 2, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
515 // SetParams( 3, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
516  }
517 //
519  if ( NumberOfTargetNucleons > 26 ) {
521  } else {
523  }
524  SetProjMinDiffMass( 1.16 ); // GeV
525  SetProjMinNonDiffMass( 1.16 ); // GeV
526  SetTarMinDiffMass( 1.16 ); // GeV
527  SetTarMinNonDiffMass( 1.16 ); // GeV
528  SetAveragePt2( 0.15 ); // GeV^2 // Uzhi Oct 2014
529  SetProbLogDistrPrD( 0.3 ); // Uzhi Oct 2014 0.5
530  SetProbLogDistr(0.3 ); // 0.5
531 
532  } else if( ProjectilePDGcode < -1000 ) { // Projectile is anti_baryon
533 
534  // Proc# A1 B1 A2 B2 A3 Atop Ymin
535  SetParams( 0, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 ); // Qexchange without Exc.
536  SetParams( 1, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 ); // Qexchange with Exc.
537  SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Projectile diffraction
538  SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Target diffraction
539  SetParams( 4, 1.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.93 ); // Qexchange with Exc. Additional multiply
540 
541  if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
542  SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
543 // SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
544  }
547  SetProjMinDiffMass( ProjectileMass + 0.22 ); // GeV
548  SetProjMinNonDiffMass( ProjectileMass + 0.22 ); // GeV
549  SetTarMinDiffMass( TargetMass + 0.22 ); // GeV
550  SetTarMinNonDiffMass( TargetMass + 0.22 ); // GeV
551  SetAveragePt2( 0.15 ); // GeV^2 // Uzhi Oct 2014
552  SetProbLogDistrPrD( 0.3 ); // Uzhi Oct 2014
553  SetProbLogDistr( 0.3 );
554 
555  } else if ( ProjectileabsPDGcode == 211 || ProjectilePDGcode == 111 ) { // Projectile is Pion
556 
557  // Proc# A1 B1 A2 B2 A3 Atop Ymin
558  SetParams( 0, 720.0, 2.5 , 2.3 , 1.0, 0., 1. , 2.7 );
559  SetParams( 1, 12.87, 0.5 , -44.91, 1.0, 0., 0. , 2.5 );
560  SetParams( 2, 0.086, 0. , -0.3 , 0.5, 0., 0. , 2.5 );
561  SetParams( 3, 32.8, 1.0 , -114.5 , 1.5, 0.084, 0. , 2.5 );
562  SetParams( 4, 1.0, 0.0 , -3.49, 0.5, 0.0, 0. , 2.5 ); // Qexchange with Exc. Additional multiply
563 
564  if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
565  SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
566 // SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
567  }
568 
569  SetDeltaProbAtQuarkExchange( 0.56 ); // (0.35)
570  SetProjMinDiffMass( 0.5 ); // (0.5) // GeV
571  SetProjMinNonDiffMass( 0.5 ); // (0.5) // GeV
572  SetTarMinDiffMass( 1.16 ); // GeV
573  SetTarMinNonDiffMass( 1.16 ); // GeV
574  SetAveragePt2( 0.15 ); // GeV^2 // Uzhi Oct 2014
575  SetProbLogDistrPrD( 0.3 ); // Uzhi Oct 2014
576  SetProbLogDistr( 0.3 );
577 
578  } else if ( ProjectileabsPDGcode == 321 || ProjectileabsPDGcode == 311 ||
579  ProjectilePDGcode == 130 || ProjectilePDGcode == 310 ) { // Projectile is Kaon
580 
581  // Proc# A1 B1 A2 B2 A3 Atop Ymin
582  SetParams( 0, 60.0 , 2.5 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Qexchange without Exc.
583  SetParams( 1, 6.0 , 1.0 , -24.33 , 2.0 , 0.0, 0.0 , 1.40 ); // Qexchange with Exc.
584  SetParams( 2, 2.76, 1.2 , -22.5 , 2.7 ,0.04, 0.0 , 1.40 ); // Projectile diffraction
585  SetParams( 3, 1.09, 0.5 , -8.88 , 2. ,0.05, 0.0 , 1.40 ); // Target diffraction
586  SetParams( 4, 1.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.93 ); // Qexchange with Exc. Additional multiply
587 
588  if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
589  SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
590 // SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
591  }
593  SetProjMinDiffMass( 0.7 ); // (1.4) // (0.7) // GeV
594  SetProjMinNonDiffMass( 0.7 ); // (1.4) // (0.7) // GeV
595  SetTarMinDiffMass( 1.16 ); // GeV
596  SetTarMinNonDiffMass( 1.16 ); // GeV
597  SetAveragePt2( 0.15 ); // GeV^2 // Uzhi Oct 2014
598  SetProbLogDistrPrD( 0.5 ); // Uzhi Oct 2014
599  SetProbLogDistr( 0.3 ); // Uzhi 5.06.2012
600 
601  } else { // Projectile is undefined, Nucleon assumed
602 
603  // Proc# A1 B1 A2 B2 A3 Atop Ymin
604  SetParams( 0, 13.71, 1.75, -214.5, 4.25, 0.0, 0.5 , 1.1 ); // Qexchange without Exc.
605  SetParams( 1, 25.0, 1.0, -50.34, 1.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc.
606  SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Projectile diffraction
607  SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Target diffraction
608  SetParams( 4, 1.0, 0.0 , -2.01 , 0.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc. Additional multiply
609 
610  if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
611  SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
612 // SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
613  }
614  SetDeltaProbAtQuarkExchange( 0.0 ); // 7 June 2011
616  SetProjMinDiffMass( ProjectileMass + 0.22 ); // GeV
617  SetProjMinNonDiffMass( ProjectileMass + 0.22 ); // GeV
618  SetTarMinDiffMass( TargetMass + 0.22 ); // GeV
619  SetTarMinNonDiffMass( TargetMass + 0.22 ); // GeV
620  SetAveragePt2( 0.15 ); // GeV^2 // Uzhi Oct 2014
621  SetProbLogDistrPrD( 0.3 ); // Uzhi Oct 2014
622  SetProbLogDistr( 0.3 );
623 
624  }
625 
626  // Set parameters of a string kink
627 // SetPt2Kink( 6.0*GeV*GeV ); // Uzhi Oct 2014
628  SetPt2Kink( 0.0*GeV*GeV ); // Uzhi Oct 2014 // Uzhi to switch off kinky strings
629  G4double Puubar( 1.0/3.0 ), Pddbar( 1.0/3.0 ), Pssbar( 1.0/3.0 ); // SU(3) symmetry
630 // G4double Puubar( 0.41 ), Pddbar( 0.41 ), Pssbar( 0.18 ); // Broken SU(3) symmetry
631  SetQuarkProbabilitiesAtGluonSplitUp( Puubar, Pddbar, Pssbar );
632 
633  // Set parameters of nuclear destruction
634  if ( ProjectileabsPDGcode < 1000 ) { // Meson projectile
635  SetMaxNumberOfCollisions( Plab, 2.0 ); // 3.0 )
636  SetCofNuclearDestruction( 1.0*std::exp( 4.0*(Ylab - 2.1) )/
637  ( 1.0 + std::exp( 4.0*(Ylab - 2.1) ) ) ); // 0.62 1.0
640  SetPt2ofNuclearDestruction( ( 0.035 + 0.04*std::exp( 4.0*(Ylab - 2.5) )/
641  ( 1.0 + std::exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV ); // 0.09
642  SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
644  } else if ( ProjectilePDGcode < -1000 ) { // for anti-baryon projectile
645  //G4cout << "Nucl destruct Anti Bar" << G4endl;
646  SetMaxNumberOfCollisions( Plab, 2.0 ); // 3.0 )
647  SetCofNuclearDestruction( 1.0*std::exp( 4.0*(Ylab - 2.1) )/
648  ( 1.0 + std::exp( 4.0*(Ylab - 2.1) ) ) ); // 0.62 1.0
651  SetPt2ofNuclearDestruction( ( 0.035 + 0.04*std::exp( 4.0*(Ylab - 2.5) )/
652  ( 1.0 + std::exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV ); // 0.09
653  SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
655  if ( Plab < 2.0 ) { // 2 GeV/c
656  // For slow anti-baryon we have to garanty putting on mass-shell
658  SetR2ofNuclearDestruction( 1.5*fermi*fermi );
659  SetDofNuclearDestruction( 0.01 );
660  SetPt2ofNuclearDestruction( 0.035*GeV*GeV );
661  SetMaxPt2ofNuclearDestruction( 0.04*GeV*GeV );
662  //SetExcitationEnergyPerWoundedNucleon( 0.0 ); // ?????
663  }
664  } else { // Projectile baryon assumed
665  SetMaxNumberOfCollisions( Plab, 2.0 ); // 3.0 )
666  SetCofNuclearDestruction( 1.0*std::exp( 4.0*(Ylab - 2.1) )/
667  ( 1.0 + std::exp( 4.0*(Ylab - 2.1) ) ) ); // 0.62 1.0
670  SetPt2ofNuclearDestruction( ( 0.035 + 0.04*std::exp( 4.0*(Ylab - 2.5) )/
671  ( 1.0 + std::exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV ); // 0.09
672  SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
674  }
675 
676  //SetCofNuclearDestruction( 0.47*std::exp( 2.0*(Ylab - 2.5) )/( 1.0 + std::exp( 2.0*(Ylab - 2.5) ) ) );
677  //SetPt2ofNuclearDestruction( ( 0.035 + 0.1*std::exp( 4.0*(Ylab - 3.0) )/( 1.0 + std::exp( 4.0*(Ylab - 3.0) ) ) )*GeV*GeV );
678 
679  //SetMagQuarkExchange( 120.0 ); // 210.0 PipP
680  //SetSlopeQuarkExchange( 2.0 );
681  //SetDeltaProbAtQuarkExchange( 0.6 );
682  //SetProjMinDiffMass( 0.7 ); // GeV 1.1
683  //SetProjMinNonDiffMass( 0.7 ); // GeV
684  //SetProbabilityOfProjDiff( 0.0); // 0.85*std::pow( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
685  //SetTarMinDiffMass( 1.1 ); // GeV
686  //SetTarMinNonDiffMass( 1.1 ); // GeV
687  //SetProbabilityOfTarDiff( 0.0 ); // 0.85*std::pow( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
688 
689 //SetAveragePt2( 0.0 ); // GeV^2 0.3
690  //------------------------------------
691 //SetProbabilityOfElasticScatt( 1.0, 1.0); //(Xtotal, Xelastic);
692  //SetProbabilityOfProjDiff( 1.0*0.62*std::pow( s/GeV/GeV, -0.51 ) ); // 0->1
693  //SetProbabilityOfTarDiff( 4.0*0.62*std::pow( s/GeV/GeV, -0.51 ) ); // 2->4
694  //SetAveragePt2( 0.3 ); // (0.15)
695 //SetAvaragePt2ofElasticScattering( 0.0 );
696 
697  //SetMaxNumberOfCollisions( Plab, 6.0 ); //(4.0*(Plab + 0.01), Plab); // 6.0 );
698  //SetAveragePt2( 0.15 );
699 // G4cout << "Cnd " << GetCofNuclearDestruction() << G4endl;
700 //SetCofNuclearDestruction( 0.0 ); // (0.2) // (0.4) 0.5
701 //SetExcitationEnergyPerWoundedNucleon( 0.0*MeV ); // (75.0*MeV)
702 //SetDofNuclearDestruction( 0.0 ); // 0.3 0.5
703 //SetPt2ofNuclearDestruction( 0.0*GeV*GeV ); // (0.168*GeV*GeV)
704  //G4cout << "Pt2 " << GetPt2ofNuclearDestruction()/GeV/GeV << G4endl;
705  //G4int Uzhi; G4cin >> Uzhi;
706 
707 }
708 
709 
710 //============================================================================
711 
713  G4double Prob( 0.0 );
714  if ( y < ProcParams[ProcN][6] ) {
715  Prob = ProcParams[ProcN][5];
716  if(Prob < 0.) Prob=0.; // Uzhi Oct 2014
717  return Prob;
718  }
719  Prob = ProcParams[ProcN][0] * std::exp( -ProcParams[ProcN][1]*y ) +
720  ProcParams[ProcN][2] * std::exp( -ProcParams[ProcN][3]*y ) +
721  ProcParams[ProcN][4];
722  if(Prob < 0.) Prob=0.; // Uzhi Oct 2014
723  return Prob;
724 }
static const double MeV
Definition: G4SIunits.hh:193
static G4ThreadLocal G4ChipsComponentXS * chipsComponentXSinstance
void SetProbLogDistr(const G4double aValue)
const G4double pi
void SetProbLogDistrPrD(const G4double aValue)
#define G4ThreadLocal
Definition: tls.hh:89
void SetTarMinNonDiffMass(const G4double aValue)
void SetDofNuclearDestruction(const G4double aValue)
void SetSlope(const G4double Slope)
G4ChipsComponentXS * FTFxsManager
int G4int
Definition: G4Types.hh:78
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
void SetGamma0(const G4double Gamma0)
G4GLOB_DLL std::ostream G4cout
void SetAveragePt2(const G4double aValue)
void SetElastisCrossSection(const G4double Xelastic)
bool G4bool
Definition: G4Types.hh:79
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()
Definition: G4Proton.cc:93
void SetProjMinNonDiffMass(const G4double aValue)
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
void SetQuarkProbabilitiesAtGluonSplitUp(const G4double Puubar, const G4double Pddbar, const G4double Pssbar)
static const double GeV
Definition: G4SIunits.hh:196
virtual G4double GetElasticElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double N)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
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()
Definition: G4PionMinus.cc:98
G4double ProcParams[5][7]
void SetCofNuclearDestruction(const G4double aValue)
static const double millibarn
Definition: G4SIunits.hh:96
void SetProbOfSameQuarkExchange(const G4double aValue)
G4double FTFXannihilation
#define G4endl
Definition: G4ios.hh:61
virtual G4double GetTotalElementCrossSection(const G4ParticleDefinition *aParticle, G4double kinEnergy, G4int Z, G4double N)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
void SetProbabilityOfAnnihilation(const G4double aValue)
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4double GetPDGCharge() const
static G4ThreadLocal bool chipsComponentXSisInitialized
void SetPt2ofNuclearDestruction(const G4double aValue)
void SetAvaragePt2ofElasticScattering(const G4double aPt2)
void SetTarMinDiffMass(const G4double aValue)
static const double fermi
Definition: G4SIunits.hh:93
void SetR2ofNuclearDestruction(const G4double aValue)