ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FTFParameters.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4FTFParameters.cc
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 //
28 
29 #include <utility>
30 
31 #include "G4FTFParameters.hh"
32 
33 #include "G4ios.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 
37 #include "G4ParticleDefinition.hh"
38 #include "G4Proton.hh"
39 #include "G4Neutron.hh"
40 #include "G4PionPlus.hh"
41 #include "G4PionMinus.hh"
42 #include "G4KaonPlus.hh"
43 #include "G4KaonMinus.hh"
44 
45 #include "G4CrossSectionDataSetRegistry.hh"
46 #include "G4VComponentCrossSection.hh"
47 #include "G4ComponentGGHadronNucleusXsc.hh"
48 #include "G4LundStringFragmentation.hh"
49 
50 #include "G4Exp.hh"
51 #include "G4Log.hh"
52 #include "G4Pow.hh"
53 
54 #include "G4HadronicDeveloperParameters.hh"
55 
56 //============================================================================
57 
58 //#define debugFTFparams
59 
60 //============================================================================
61 
62 G4FTFParameters::G4FTFParameters()
63 {
64  StringMass = new G4LundStringFragmentation; // for estimation of min. mass of diffr. states
65  Reset();
66  csGGinstance =
67  G4CrossSectionDataSetRegistry::Instance()->GetComponentCrossSection("Glauber-Gribov");
68  if (!csGGinstance) {
69  csGGinstance = new G4ComponentGGHadronNucleusXsc();
70  }
71 
72  // Set parameters of a string kink
73  SetPt2Kink( 0.0*GeV*GeV ); // To switch off kinky strings (bad results obtained with 6.0*GeV*GeV)
74  G4double Puubar( 1.0/3.0 ), Pddbar( 1.0/3.0 ), Pssbar( 1.0/3.0 ); // SU(3) symmetry
75  //G4double Puubar( 0.41 ), Pddbar( 0.41 ), Pssbar( 0.18 ); // Broken SU(3) symmetry
76  SetQuarkProbabilitiesAtGluonSplitUp( Puubar, Pddbar, Pssbar );
77 }
78 
79 //============================================================================
80 
81 void G4FTFParameters::InitForInteraction( const G4ParticleDefinition* particle,
82  G4int theA, G4int theZ, G4double PlabPerParticle )
83 {
84  Reset();
85 
86  G4int ProjectilePDGcode = particle->GetPDGEncoding();
87  G4int ProjectileabsPDGcode = std::abs( ProjectilePDGcode );
88  G4double ProjectileMass = particle->GetPDGMass();
89  G4double ProjectileMass2 = ProjectileMass * ProjectileMass;
90 
91  G4int ProjectileBaryonNumber( 0 ), AbsProjectileBaryonNumber( 0 ), AbsProjectileCharge( 0 );
92  G4bool ProjectileIsNucleus = false;
93 
94  if ( std::abs( particle->GetBaryonNumber() ) > 1 ) { // The projectile is a nucleus
95  ProjectileIsNucleus = true;
96  ProjectileBaryonNumber = particle->GetBaryonNumber();
97  AbsProjectileBaryonNumber = std::abs( ProjectileBaryonNumber );
98  AbsProjectileCharge = std::abs( G4int( particle->GetPDGCharge() ) );
99  if ( ProjectileBaryonNumber > 1 ) {
100  ProjectilePDGcode = 2212; ProjectileabsPDGcode = 2212; // Proton
101  } else {
102  ProjectilePDGcode = -2212; ProjectileabsPDGcode = 2212; // Anti-Proton
103  }
104  ProjectileMass = G4Proton::Proton()->GetPDGMass();
105  ProjectileMass2 = sqr( ProjectileMass );
106  }
107 
108  G4double TargetMass = G4Proton::Proton()->GetPDGMass();
109  G4double TargetMass2 = TargetMass * TargetMass;
110 
111  G4double Plab = PlabPerParticle;
112  G4double Elab = std::sqrt( Plab*Plab + ProjectileMass2 );
113  G4double KineticEnergy = Elab - ProjectileMass;
114 
115  G4double S = ProjectileMass2 + TargetMass2 + 2.0*TargetMass*Elab;
116 
117  #ifdef debugFTFparams
118  G4cout << "--------- FTF Parameters --------------" << G4endl << "Proj Plab "
119  << ProjectilePDGcode << " " << Plab << G4endl << "Mass KinE " << ProjectileMass
120  << " " << KineticEnergy << G4endl << " A Z " << theA << " " << theZ << G4endl;
121  #endif
122 
123  G4double Ylab, Xtotal( 0.0 ), Xelastic( 0.0 ), Xannihilation( 0.0 );
124  G4int NumberOfTargetNucleons;
125 
126  Ylab = 0.5 * G4Log( (Elab + Plab)/(Elab - Plab) );
127 
128  G4double ECMSsqr = S/GeV/GeV;
129  G4double SqrtS = std::sqrt( S )/GeV;
130 
131  #ifdef debugFTFparams
132  G4cout << "Sqrt(s) " << SqrtS << G4endl;
133  #endif
134 
135  TargetMass /= GeV; TargetMass2 /= (GeV*GeV);
136  ProjectileMass /= GeV; ProjectileMass2 /= (GeV*GeV);
137 
138  Plab /= GeV;
139  G4double Xftf = 0.0;
140 
141  G4int NumberOfTargetProtons = theZ;
142  G4int NumberOfTargetNeutrons = theA - theZ;
143  NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
144 
145  // ---------- hadron projectile ----------------
146  if ( AbsProjectileBaryonNumber <= 1 ) { // Projectile is hadron or baryon
147 
148  // Interaction on P
149  G4double xTtP = csGGinstance->GetTotalIsotopeCrossSection( particle, KineticEnergy, 1, 1);
150  G4double xElP = csGGinstance->GetElasticIsotopeCrossSection(particle, KineticEnergy, 1, 1);
151 
152  // Interaction on N
153  G4double xTtN = csGGinstance->GetTotalIsotopeCrossSection( particle, KineticEnergy, 0, 1);
154  G4double xElN = csGGinstance->GetElasticIsotopeCrossSection(particle, KineticEnergy, 0, 1);
155 
156  // Average properties of h+N interactions
157  Xtotal = ( NumberOfTargetProtons * xTtP + NumberOfTargetNeutrons * xTtN ) / NumberOfTargetNucleons;
158  Xelastic = ( NumberOfTargetProtons * xElP + NumberOfTargetNeutrons * xElN ) / NumberOfTargetNucleons;
159  Xannihilation = 0.0;
160 
161  Xtotal /= millibarn;
162  Xelastic /= millibarn;
163 
164  #ifdef debugFTFparams
165  G4cout<<"Estimated cross sections (total and elastic) of h+N interactions "<<Xtotal<<" "<<Xelastic<<" (mb)"<<G4endl;
166  #endif
167  }
168 
169  // ---------- nucleus projectile ----------------
170  if ( ProjectileIsNucleus && ProjectileBaryonNumber > 1 ) {
171 
172  #ifdef debugFTFparams
173  G4cout<<"Projectile is a nucleus: A and Z - "<<ProjectileBaryonNumber<<" "<<ProjectileCharge<<G4endl;
174  #endif
175 
176  const G4ParticleDefinition* Proton = G4Proton::Proton();
177  // Interaction on P
178  G4double XtotPP = csGGinstance->GetTotalIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
179  G4double XelPP = csGGinstance->GetElasticIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
180 
181  const G4ParticleDefinition* Neutron = G4Neutron::Neutron();
182  // Interaction on N
183  G4double XtotPN = csGGinstance->GetTotalIsotopeCrossSection(Neutron, KineticEnergy, 0, 1);
184  G4double XelPN = csGGinstance->GetElasticIsotopeCrossSection(Neutron, KineticEnergy, 0, 1);
185 
186  #ifdef debugFTFparams
187  G4cout << "XsPP (total and elastic) " << XtotPP/millibarn << " " << XelPP/millibarn <<" (mb)"<< G4endl
188  << "XsPN (total and elastic) " << XtotPN/millibarn << " " << XelPN/millibarn <<" (mb)"<< G4endl;
189  #endif
190 
191  Xtotal = (
192  AbsProjectileCharge * NumberOfTargetProtons * XtotPP +
193  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
194  NumberOfTargetNeutrons * XtotPP
195  +
196  ( AbsProjectileCharge * NumberOfTargetNeutrons +
197  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
198  NumberOfTargetProtons ) * XtotPN
199  ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
200  Xelastic= (
201  AbsProjectileCharge * NumberOfTargetProtons * XelPP +
202  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
203  NumberOfTargetNeutrons * XelPP
204  +
205  ( AbsProjectileCharge * NumberOfTargetNeutrons +
206  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
207  NumberOfTargetProtons ) * XelPN
208  ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
209 
210  Xannihilation = 0.0;
211  Xtotal /= millibarn;
212  Xelastic /= millibarn;
213  }
214 
215  // ---------- The projectile is anti-baryon or anti-nucleus ----------------
216  // anti Sigma^0_c anti Delta^-
217  if ( ProjectilePDGcode >= -4112 && ProjectilePDGcode <= -1114 ) {
218  // Only non-strange and strange baryons are considered
219 
220  #ifdef debugFTFparams
221  G4cout<<"Projectile is a anti-baryon or anti-nucleus - "<<ProjectileBaryonNumber<<" "<<ProjectileCharge<<G4endl;
222  G4cout<<"(Only non-strange and strange baryons are considered)"<<G4endl;
223  #endif
224 
225  G4double X_a( 0.0 ), X_b( 0.0 ), X_c( 0.0 ), X_d( 0.0 );
226  G4double MesonProdThreshold = ProjectileMass + TargetMass +
227  ( 2.0 * 0.14 + 0.016 ); // 2 Mpi + DeltaE;
228 
229  if ( PlabPerParticle < 40.0*MeV ) { // Low energy limits. Projectile at rest.
230  Xtotal = 1512.9; // mb
231  Xelastic = 473.2; // mb
232  X_a = 625.1; // mb
233  X_b = 9.780; // mb
234  X_c = 49.989; // mb
235  X_d = 6.614; // mb
236  } else { // Total and elastic cross section of PbarP interactions a'la Arkhipov
237  G4double LogS = G4Log( ECMSsqr / 33.0625 );
238  G4double Xasmpt = 36.04 + 0.304*LogS*LogS; // mb
239  LogS = G4Log( SqrtS / 20.74 );
240  G4double Basmpt = 11.92 + 0.3036*LogS*LogS; // GeV^(-2)
241  G4double R0 = std::sqrt( 0.40874044*Xasmpt - Basmpt ); // GeV^(-1)
242 
243  G4double FlowF = SqrtS / std::sqrt( ECMSsqr*ECMSsqr + ProjectileMass2*ProjectileMass2 +
244  TargetMass2*TargetMass2 - 2.0*ECMSsqr*ProjectileMass2
245  - 2.0*ECMSsqr*TargetMass2
246  - 2.0*ProjectileMass2*TargetMass2 );
247 
248  Xtotal = Xasmpt * ( 1.0 + 13.55*FlowF/R0/R0/R0*
249  (1.0 - 4.47/SqrtS + 12.38/ECMSsqr - 12.43/SqrtS/ECMSsqr) ); // mb
250 
251  Xasmpt = 4.4 + 0.101*LogS*LogS; // mb
252  Xelastic = Xasmpt * ( 1.0 + 59.27*FlowF/R0/R0/R0*
253  (1.0 - 6.95/SqrtS + 23.54/ECMSsqr - 25.34/SqrtS/ECMSsqr ) ); // mb
254 
255  X_a = 25.0*FlowF; // mb, 3-shirts diagram
256 
257  if ( SqrtS < MesonProdThreshold ) {
258  X_b = 3.13 + 140.0*G4Pow::GetInstance()->powA( MesonProdThreshold - SqrtS, 2.5 ); // mb anti-quark-quark annihilation
259  Xelastic -= 3.0*X_b; // Xel-X(PbarP->NNbar)
260  } else {
261  X_b = 6.8/SqrtS; // mb anti-quark-quark annihilation
262  Xelastic -= 3.0*X_b; // Xel-X(PbarP->NNbar)
263  }
264 
265  X_c = 2.0*FlowF*sqr( ProjectileMass + TargetMass )/ECMSsqr; // mb rearrangement
266 
267  X_d = 23.3/ECMSsqr; // mb anti-quark-quark string creation
268  }
269 
270  G4double Xann_on_P( 0.0), Xann_on_N( 0.0 );
271 
272  if ( ProjectilePDGcode == -2212 ) { // Pbar+P/N
273  Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0;
274  Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
275  } else if ( ProjectilePDGcode == -2112 ) { // NeutrBar+P/N
276  Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
277  Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0;
278  } else if ( ProjectilePDGcode == -3122 ) { // LambdaBar+P/N
279  Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
280  Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
281  } else if ( ProjectilePDGcode == -3112 ) { // Sigma-Bar+P/N
282  Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
283  Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
284  } else if ( ProjectilePDGcode == -3212 ) { // Sigma0Bar+P/N
285  Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
286  Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
287  } else if ( ProjectilePDGcode == -3222 ) { // Sigma+Bar+P/N
288  Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
289  Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
290  } else if ( ProjectilePDGcode == -3312 ) { // Xi-Bar+P/N
291  Xann_on_P = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
292  Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
293  } else if ( ProjectilePDGcode == -3322 ) { // Xi0Bar+P/N
294  Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
295  Xann_on_N = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
296  } else if ( ProjectilePDGcode == -3334 ) { // Omega-Bar+P/N
297  Xann_on_P = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
298  Xann_on_N = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
299  } else {
300  G4cout << "Unknown anti-baryon for FTF annihilation" << G4endl;
301  }
302 
303  //G4cout << "Sum " << Xann_on_P << G4endl;
304 
305  if ( ! ProjectileIsNucleus ) { // Projectile is anti-baryon
306  Xannihilation = ( NumberOfTargetProtons * Xann_on_P + NumberOfTargetNeutrons * Xann_on_N )
307  / NumberOfTargetNucleons;
308  } else { // Projectile is a nucleus
309  Xannihilation = (
310  ( AbsProjectileCharge * NumberOfTargetProtons +
311  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
312  NumberOfTargetNeutrons ) * Xann_on_P
313  +
314  ( AbsProjectileCharge * NumberOfTargetNeutrons +
315  ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
316  NumberOfTargetProtons ) * Xann_on_N
317  ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
318  }
319 
320  //G4double Xftf = 0.0;
321  MesonProdThreshold = ProjectileMass + TargetMass + (0.14 + 0.08); // Mpi + DeltaE
322  if ( SqrtS > MesonProdThreshold ) {
323  Xftf = 36.0 * ( 1.0 - MesonProdThreshold/SqrtS );
324  }
325 
326  Xtotal = Xelastic + Xannihilation + Xftf;
327 
328  #ifdef debugFTFparams
329  G4cout << "Plab Xtotal, Xelastic Xinel Xftf " << Plab << " " << Xtotal << " " << Xelastic
330  << " " << Xtotal - Xelastic << " " << Xtotal - Xelastic - Xannihilation << " (mb)"<< G4endl
331  << "Plab Xelastic/Xtotal, Xann/Xin " << Plab << " " << Xelastic/Xtotal << " "
332  << Xannihilation/(Xtotal - Xelastic) << G4endl;
333  #endif
334 
335  }
336 
337  if ( Xtotal == 0.0 ) { // Projectile is undefined, Nucleon assumed
338 
339  const G4ParticleDefinition* Proton = G4Proton::Proton();
340  // Interaction on P
341  G4double XtotPP = csGGinstance->GetTotalIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
342  G4double XelPP = csGGinstance->GetElasticIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
343 
344  // Interaction on N
345  G4double XtotPN = csGGinstance->GetTotalIsotopeCrossSection(Proton, KineticEnergy, 0, 1);
346  G4double XelPN = csGGinstance->GetElasticIsotopeCrossSection(Proton, KineticEnergy, 0, 1);
347 
348  Xtotal = ( NumberOfTargetProtons * XtotPP + NumberOfTargetNeutrons * XtotPN )
349  / NumberOfTargetNucleons;
350  Xelastic = ( NumberOfTargetProtons * XelPP + NumberOfTargetNeutrons * XelPN )
351  / NumberOfTargetNucleons;
352  Xannihilation = 0.0;
353  Xtotal /= millibarn;
354  Xelastic /= millibarn;
355  };
356 
357  // Geometrical parameters
358  SetTotalCrossSection( Xtotal );
359  SetElastisCrossSection( Xelastic );
360  SetInelasticCrossSection( Xtotal - Xelastic );
361 
362  // Interactions with elastic and inelastic collisions
363  SetProbabilityOfElasticScatt( Xtotal, Xelastic );
364 
365  SetRadiusOfHNinteractions2( Xtotal/pi/10.0 );
366 
367  if ( ( Xtotal - Xelastic ) == 0.0 ) {
368  SetProbabilityOfAnnihilation( 0.0 );
369  } else {
370  SetProbabilityOfAnnihilation( Xannihilation / (Xtotal - Xelastic) );
371  }
372 
373  if(Xelastic > 0.0) {
374  SetSlope( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 );// Slope parameter of elastic scattering
375  // (GeV/c)^(-2))
376  // Parameters of elastic scattering
377  // Gaussian parametrization of elastic scattering amplitude assumed
378  SetAvaragePt2ofElasticScattering( 1.0/( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 )*GeV*GeV );
379  } else {
380  SetSlope(1.0);
381  SetAvaragePt2ofElasticScattering( 0.0);
382  }
383  SetGamma0( GetSlope()*Xtotal/10.0/2.0/pi );
384 
385  G4double Xinel = Xtotal - Xelastic;
386 
387  #ifdef debugFTFparams
388  G4cout<< "Slope of hN elastic scattering" << GetSlope() << G4endl;
389  G4cout << "AvaragePt2ofElasticScattering " << GetAvaragePt2ofElasticScattering() << G4endl;
390  G4cout<<"Parameters of excitation for projectile "<<ProjectilePDGcode<< G4endl;
391  #endif
392 
393  if ( (ProjectilePDGcode == 2212) || (ProjectilePDGcode == 2112) ) { // Projectile is proton or neutron
394 
395  // A process probability is parameterized as Prob = A_1*exp(-A_2*y) + A_3*exp(-A_4*y) + A_top
396  // y is a rapidity of a partcle in the target nucleus. Ymin is a minimal rapidity below it X=0
397 
398  // Proc# A1 B1 A2 B2 A3 Atop Ymin
399  /* original hadr-string-diff-V10-03-07 (similar to 10.3.x)
400  SetParams( 0, 13.71, 1.75, -214.5, 4.25, 0.0, 0.5 , 1.1 ); // Qexchange without Exc.
401  SetParams( 1, 25.0, 1.0, -50.34, 1.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc.
402  SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Projectile diffraction
403  SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Target diffraction
404  SetParams( 4, 1.0, 0.0 , -2.01 , 0.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc. Additional multiplier
405  */
406 
407  // Proc#
408  SetParams( 0, fParCollBaryonProj.GetProc0A1(), fParCollBaryonProj.GetProc0B1(),
409  fParCollBaryonProj.GetProc0A2(), fParCollBaryonProj.GetProc0B2(),
410  fParCollBaryonProj.GetProc0A3(), fParCollBaryonProj.GetProc0Atop(),
411  fParCollBaryonProj.GetProc0Ymin() ); // Qexchange without Exc.
412  SetParams( 1, fParCollBaryonProj.GetProc1A1(), fParCollBaryonProj.GetProc1B1(),
413  fParCollBaryonProj.GetProc1A2(), fParCollBaryonProj.GetProc1B2(),
414  fParCollBaryonProj.GetProc1A3(), fParCollBaryonProj.GetProc1Atop(),
415  fParCollBaryonProj.GetProc1Ymin() ); // Qexchange with Exc.
416  if ( Xinel > 0.0 ) {
417  SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 ); // Projectile diffraction
418  SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 ); // Target diffraction
419 
420  SetParams( 4, fParCollBaryonProj.GetProc4A1(), fParCollBaryonProj.GetProc4B1(),
421  fParCollBaryonProj.GetProc4A2(), fParCollBaryonProj.GetProc4B2(),
422  fParCollBaryonProj.GetProc4A3(), fParCollBaryonProj.GetProc4Atop(),
423  fParCollBaryonProj.GetProc4Ymin() ); // Qexchange with Exc. Additional multiplier
424  } else { // if Xinel=0., zero everything out (obviously)
425  SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
426  SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
427  SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
428  }
429 
430  if ( AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10 ) {
431  //
432  // It is not decided what to do with diffraction dissociation in Had-Nucl and Nucl-Nucl interactions
433  // For the moment both ProjDiffDisso & TgtDiffDisso for A > 10 are set to false,
434  // so both projectile and target diffraction are turned OFF
435  //
436  if ( ! fParCollBaryonProj.IsProjDiffDissociation() )
437  SetParams( 2, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
438  if ( ! fParCollBaryonProj.IsTgtDiffDissociation() )
439  SetParams( 3, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
440  }
441 
442  SetDeltaProbAtQuarkExchange( fParCollBaryonProj.GetDeltaProbAtQuarkExchange() );
443 
444  if ( NumberOfTargetNucleons > 26 ) {
445  SetProbOfSameQuarkExchange( 1.0 );
446  } else {
447  SetProbOfSameQuarkExchange( fParCollBaryonProj.GetProbOfSameQuarkExchange() );
448  }
449 
450  SetProjMinDiffMass( fParCollBaryonProj.GetProjMinDiffMass() ); // GeV
451  SetProjMinNonDiffMass( fParCollBaryonProj.GetProjMinNonDiffMass() ); // GeV
452 
453  SetTarMinDiffMass( fParCollBaryonProj.GetTgtMinDiffMass() ); // GeV
454  SetTarMinNonDiffMass( fParCollBaryonProj.GetTgtMinNonDiffMass() ); // GeV
455 
456  SetAveragePt2( fParCollBaryonProj.GetAveragePt2() ); // GeV^2
457  SetProbLogDistrPrD( fParCollBaryonProj.GetProbLogDistrPrD() );
458  SetProbLogDistr( fParCollBaryonProj.GetProbLogDistr() );
459 
460  } else if ( ProjectilePDGcode == -2212 || ProjectilePDGcode == -2112 ) { // Projectile is anti_proton or anti_neutron
461 
462  // Proc# A1 B1 A2 B2 A3 Atop Ymin
463  SetParams( 0, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 ); // Qexchange without Exc.
464  SetParams( 1, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 ); // Qexchange with Exc.
465  if ( Xinel > 0.) {
466  SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 ); // Projectile diffraction
467  SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 ); // Target diffraction
468  SetParams( 4, 1.0, 0.0 , 0.0, 0.0 , 0.0, 0.0 , 0.93 ); // Qexchange with Exc. Additional multiply
469  } else {
470  SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0 );
471  SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0 );
472  SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0 );
473  }
474 
475  if ( AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10 ) {
476  //
477  // It is not decided what to do with diffraction dissociation in Had-Nucl and Nucl-Nucl interactions
478  // For the moment both ProjDiffDisso & TgtDiffDisso are set to false,
479  // so both projectile and target diffraction are turned OFF
480  //
481  if ( ! fParCollBaryonProj.IsProjDiffDissociation() )
482  SetParams( 2, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
483  if ( ! fParCollBaryonProj.IsTgtDiffDissociation() )
484  SetParams( 3, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
485  }
486 
487  SetDeltaProbAtQuarkExchange( 0.0 );
488  SetProbOfSameQuarkExchange( 0.0 );
489  SetProjMinDiffMass( ProjectileMass + 0.22 ); // GeV
490  SetProjMinNonDiffMass( ProjectileMass + 0.22 ); // GeV
491  SetTarMinDiffMass( TargetMass + 0.22 ); // GeV
492  SetTarMinNonDiffMass( TargetMass + 0.22 ); // GeV
493  SetAveragePt2( 0.3 ); // GeV^2
494  SetProbLogDistrPrD( 0.55 );
495  SetProbLogDistr( 0.55 );
496 
497  } else if ( ProjectileabsPDGcode == 211 || ProjectilePDGcode == 111 ) { // Projectile is Pion
498 
499  // Proc# A1 B1 A2 B2 A3 Atop Ymin
500  /* --> original code
501  SetParams( 0, 150.0, 1.8 , -247.3, 2.3, 0., 1. , 2.3 );
502  SetParams( 1, 5.77, 0.6 , -5.77, 0.8, 0., 0. , 0.0 );
503  SetParams( 2, 2.27, 0.5 , -98052.0, 4.0, 0., 0. , 3.0 );
504  SetParams( 3, 7.0, 0.9, -85.28, 1.9, 0.08, 0. , 2.2 );
505  SetParams( 4, 1.0, 0.0 , -11.02, 1.0, 0.0, 0. , 2.4 ); // Qexchange with Exc. Additional multiply
506  */
507  // Proc#
508  SetParams( 0, fParCollPionProj.GetProc0A1(), fParCollPionProj.GetProc0B1(),
509  fParCollPionProj.GetProc0A2(), fParCollPionProj.GetProc0B2(),
510  fParCollPionProj.GetProc0A3(), fParCollPionProj.GetProc0Atop(),
511  fParCollPionProj.GetProc0Ymin() ); // Qexchange without Exc.
512  SetParams( 1, fParCollPionProj.GetProc1A1(), fParCollPionProj.GetProc1B1(),
513  fParCollPionProj.GetProc1A2(), fParCollPionProj.GetProc1B2(),
514  fParCollPionProj.GetProc1A3(), fParCollPionProj.GetProc1Atop(),
515  fParCollPionProj.GetProc1Ymin() ); // Qexchange with Exc.
516  SetParams( 2, fParCollPionProj.GetProc2A1(), fParCollPionProj.GetProc2B1(),
517  fParCollPionProj.GetProc2A2(), fParCollPionProj.GetProc2B2(),
518  fParCollPionProj.GetProc2A3(), fParCollPionProj.GetProc2Atop(),
519  fParCollPionProj.GetProc2Ymin() ); // Projectile diffraction
520  SetParams( 3, fParCollPionProj.GetProc3A1(), fParCollPionProj.GetProc3B1(),
521  fParCollPionProj.GetProc3A2(), fParCollPionProj.GetProc3B2(),
522  fParCollPionProj.GetProc3A3(), fParCollPionProj.GetProc3Atop(),
523  fParCollPionProj.GetProc3Ymin() ); // Target diffraction
524  SetParams( 4, fParCollPionProj.GetProc4A1(), fParCollPionProj.GetProc4B1(),
525  fParCollPionProj.GetProc4A2(), fParCollPionProj.GetProc4B2(),
526  fParCollPionProj.GetProc4A3(), fParCollPionProj.GetProc4Atop(),
527  fParCollPionProj.GetProc4Ymin() ); // Qexchange with Exc. Additional multiply
528 
529  // NOTE: how can it be |ProjectileBaryonNumber| > 10 if projectile is a pion ???
530  //
531  if ( AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10 ) {
532  if ( ! fParCollPionProj.IsProjDiffDissociation() )
533  SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
534  if ( ! fParCollPionProj.IsTgtDiffDissociation() )
535  SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
536  }
537 
538  /* original code -->
539  SetDeltaProbAtQuarkExchange( 0.56 );
540  SetProjMinDiffMass( 1.0 ); // GeV
541  SetProjMinNonDiffMass( 1.0 ); // GeV
542  SetTarMinDiffMass( 1.16 ); // GeV
543  SetTarMinNonDiffMass( 1.16 ); // GeV
544  SetAveragePt2( 0.3 ); // GeV^2
545  SetProbLogDistrPrD( 0.55 );
546  SetProbLogDistr( 0.55 );
547  */
548 
549  // JVY update, Aug.8, 2018 --> Feb.14, 2019
550  //
551  SetDeltaProbAtQuarkExchange( fParCollPionProj.GetDeltaProbAtQuarkExchange() );
552  SetProjMinDiffMass( fParCollPionProj.GetProjMinDiffMass() ); // GeV
553  SetProjMinNonDiffMass( fParCollPionProj.GetProjMinNonDiffMass() ); // GeV
554  SetTarMinDiffMass( fParCollPionProj.GetTgtMinDiffMass() ); // GeV
555  SetTarMinNonDiffMass( fParCollPionProj.GetTgtMinNonDiffMass() ); // GeV
556  SetAveragePt2( fParCollPionProj.GetAveragePt2() ); // GeV^2
557  SetProbLogDistrPrD( fParCollPionProj.GetProbLogDistrPrD() );
558  SetProbLogDistr( fParCollPionProj.GetProbLogDistr() );
559 
560  // ---> end update
561 
562  } else if ( ProjectileabsPDGcode == 321 || ProjectileabsPDGcode == 311 ||
563  ProjectilePDGcode == 130 || ProjectilePDGcode == 310 ) { // Projectile is Kaon
564 
565  // Proc# A1 B1 A2 B2 A3 Atop Ymin
566  SetParams( 0, 60.0 , 2.5 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Qexchange without Exc.
567  SetParams( 1, 6.0 , 1.0 , -24.33 , 2.0 , 0.0, 0.0 , 1.40 ); // Qexchange with Exc.
568  SetParams( 2, 2.76, 1.2 , -22.5 , 2.7 ,0.04, 0.0 , 1.40 ); // Projectile diffraction
569  SetParams( 3, 1.09, 0.5 , -8.88 , 2. ,0.05, 0.0 , 1.40 ); // Target diffraction
570  SetParams( 4, 1.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.93 ); // Qexchange with Exc. Additional multiply
571  if ( AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10 ) {
572  SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
573  SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
574  }
575 
576  SetDeltaProbAtQuarkExchange( 0.6 );
577  SetProjMinDiffMass( 0.7 ); // GeV
578  SetProjMinNonDiffMass( 0.7 ); // GeV
579  SetTarMinDiffMass( 1.16 ); // GeV
580  SetTarMinNonDiffMass( 1.16 ); // GeV
581  SetAveragePt2( 0.3 ); // GeV^2
582  SetProbLogDistrPrD( 0.55 );
583  SetProbLogDistr( 0.55 );
584 
585  } else { // Projectile is not p, n, Pi0, Pi+, Pi-, K+, K-, K0 or their anti-particles
586 
587  if ( ProjectileabsPDGcode > 1000 ) { // The projectile is a baryon as P or N
588  // Proc# A1 B1 A2 B2 A3 Atop Ymin
589  SetParams( 0, 13.71, 1.75, -30.69, 3.0 , 0.0, 1.0 , 0.93 ); // Qexchange without Exc.
590  SetParams( 1, 25.0, 1.0, -50.34, 1.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc.
591  if ( Xinel > 0.) {
592  SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Projectile diffraction
593  SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Target diffraction
594  SetParams( 4, 1.0, 0.0 , -2.01 , 0.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc. Additional multiply
595  } else {
596  SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
597  SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
598  SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
599  }
600 
601  } else { // The projectile is a meson as K+-0
602  // Proc# A1 B1 A2 B2 A3 Atop Ymin
603  SetParams( 0, 60.0 , 2.5 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Qexchange without Exc.
604  SetParams( 1, 6.0 , 1.0 , -24.33 , 2.0 , 0.0, 0.0 , 1.40 ); // Qexchange with Exc.
605  SetParams( 2, 2.76, 1.2 , -22.5 , 2.7 ,0.04, 0.0 , 1.40 ); // Projectile diffraction
606  SetParams( 3, 1.09, 0.5 , -8.88 , 2. ,0.05, 0.0 , 1.40 ); // Target diffraction
607  SetParams( 4, 1.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.93 ); // Qexchange with Exc. Additional multiply
608  }
609 
610  if ( AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10 ) {
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 
615  SetDeltaProbAtQuarkExchange( 0.0 );
616  SetProbOfSameQuarkExchange( 0.0 );
617 
618  SetProjMinDiffMass( GetMinMass(particle)/GeV );
619  SetProjMinNonDiffMass( GetMinMass(particle)/GeV );
620 
621  const G4ParticleDefinition* Neutron = G4Neutron::Neutron();
622  SetTarMinDiffMass( GetMinMass(Neutron)/GeV );
623  SetTarMinNonDiffMass( GetMinMass(Neutron)/GeV );
624 
625  SetAveragePt2( 0.3 ); // GeV^2
626  SetProbLogDistrPrD( 0.55 );
627  SetProbLogDistr( 0.55 );
628 
629  }
630 
631  #ifdef debugFTFparams
632  G4cout<<"DeltaProbAtQuarkExchange "<< GetDeltaProbAtQuarkExchange() << G4endl;
633  G4cout<<"ProbOfSameQuarkExchange "<< GetProbOfSameQuarkExchange() << G4endl;
634  G4cout<<"ProjMinDiffMass "<< GetProjMinDiffMass()/GeV <<" GeV"<< G4endl;
635  G4cout<<"ProjMinNonDiffMass "<< GetProjMinNonDiffMass() <<" GeV"<< G4endl;
636  G4cout<<"TarMinDiffMass "<< GetTarMinDiffMass() <<" GeV"<< G4endl;
637  G4cout<<"TarMinNonDiffMass "<< GetTarMinNonDiffMass() <<" GeV"<< G4endl;
638  G4cout<<"AveragePt2 "<< GetAveragePt2() <<" GeV^2"<< G4endl;
639  G4cout<<"ProbLogDistrPrD "<< GetProbLogDistrPrD() << G4endl;
640  G4cout<<"ProbLogDistrTrD "<< GetProbLogDistr() << G4endl;
641  #endif
642 
643  // Set parameters of nuclear destruction
644 
645  if ( ProjectileabsPDGcode < 1000 ) { // Meson projectile
646 
647  SetMaxNumberOfCollisions( Plab, 2.0 ); // 3.0 )
648  //
649  // target destruction
650  //
651  /* original code --->
652  SetCofNuclearDestruction( 0.00481*G4double(NumberOfTargetNucleons)*
653  G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
654 
655  SetR2ofNuclearDestruction( 1.5*fermi*fermi );
656  SetDofNuclearDestruction( 0.3 );
657  SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
658  ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
659  SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
660  SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
661  */
662  double coeff = fParCollMesonProj.GetNuclearTgtDestructP1();
663  //
664  // NOTE (JVY): Set this switch to false/true on line 138
665  //
666  if ( fParCollMesonProj.IsNuclearTgtDestructP1_ADEP() )
667  {
668  coeff *= G4double(NumberOfTargetNucleons);
669  }
670  double exfactor = G4Exp( fParCollMesonProj.GetNuclearTgtDestructP2()
671  * (Ylab-fParCollMesonProj.GetNuclearTgtDestructP3()) );
672  coeff *= exfactor;
673  coeff /= ( 1.+ exfactor );
674 
675  SetCofNuclearDestruction( coeff );
676 
677  SetR2ofNuclearDestruction( fParCollMesonProj.GetR2ofNuclearDestruct() );
678  SetDofNuclearDestruction( fParCollMesonProj.GetDofNuclearDestruct() );
679  coeff = fParCollMesonProj.GetPt2NuclearDestructP2();
680  exfactor = G4Exp( fParCollMesonProj.GetPt2NuclearDestructP3()
681  * (Ylab-fParCollMesonProj.GetPt2NuclearDestructP4()) );
682  coeff *= exfactor;
683  coeff /= ( 1. + exfactor );
684  SetPt2ofNuclearDestruction( (fParCollMesonProj.GetPt2NuclearDestructP1()+coeff)*CLHEP::GeV*CLHEP::GeV );
685 
686  SetMaxPt2ofNuclearDestruction( fParCollMesonProj.GetMaxPt2ofNuclearDestruct() );
687  SetExcitationEnergyPerWoundedNucleon( fParCollMesonProj.GetExciEnergyPerWoundedNucleon() );
688 
689  } else if ( ProjectilePDGcode == -2212 || ProjectilePDGcode == -2112 ) { // for anti-baryon projectile
690 
691  SetMaxNumberOfCollisions( Plab, 2.0 );
692 
693  SetCofNuclearDestruction( 0.00481*G4double(NumberOfTargetNucleons)*
694  G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
695  SetR2ofNuclearDestruction( 1.5*fermi*fermi );
696  SetDofNuclearDestruction( 0.3 );
697  SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
698  ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
699  SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
700  SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
701  if ( Plab < 2.0 ) { // 2 GeV/c
702  // For slow anti-baryon we have to garanty putting on mass-shell
703  SetCofNuclearDestruction( 0.0 );
704  SetR2ofNuclearDestruction( 1.5*fermi*fermi ); // this is equivalent to setting a few line above
705  // is it even necessary ?
706  SetDofNuclearDestruction( 0.01 );
707  SetPt2ofNuclearDestruction( 0.035*GeV*GeV );
708  SetMaxPt2ofNuclearDestruction( 0.04*GeV*GeV );
709  }
710 
711  } else { // Projectile baryon assumed
712 
713  // NOTE (JVY) FIXME !!! Will decide later how/if to make this one configurable...
714  //
715  SetMaxNumberOfCollisions( Plab, 2.0 );
716 
717  // projectile destruction - does NOT really matter for particle projectile, only for a nucleus projectile
718  //
719  double coeff = 0.;
720  coeff = fParCollBaryonProj.GetNuclearProjDestructP1();
721  //
722  // NOTE (JVY): Set this switch to false/true on line 136
723  //
724  if ( fParCollBaryonProj.IsNuclearProjDestructP1_NBRNDEP() )
725  {
726  coeff *= G4double(AbsProjectileBaryonNumber);
727  }
728  double exfactor = G4Exp( fParCollBaryonProj.GetNuclearProjDestructP2()*(Ylab-fParCollBaryonProj.GetNuclearProjDestructP3()) );
729  coeff *= exfactor;
730  coeff /= ( 1.+ exfactor );
731  SetCofNuclearDestructionPr( coeff );
732 
733  // target desctruction
734  //
735  coeff = fParCollBaryonProj.GetNuclearTgtDestructP1();
736  //
737  // NOTE (JVY): Set this switch to false/true on line 138
738  //
739  if ( fParCollBaryonProj.IsNuclearTgtDestructP1_ADEP() )
740  {
741  coeff *= G4double(NumberOfTargetNucleons);
742  }
743  exfactor = G4Exp( fParCollBaryonProj.GetNuclearTgtDestructP2()*(Ylab-fParCollBaryonProj.GetNuclearTgtDestructP3()) );
744  coeff *= exfactor;
745  coeff /= ( 1.+ exfactor );
746  SetCofNuclearDestruction( coeff );
747 
748  SetR2ofNuclearDestruction( fParCollBaryonProj.GetR2ofNuclearDestruct() );
749  SetDofNuclearDestruction( fParCollBaryonProj.GetDofNuclearDestruct() );
750 
751  coeff = fParCollBaryonProj.GetPt2NuclearDestructP2();
752  exfactor = G4Exp( fParCollBaryonProj.GetPt2NuclearDestructP3()*(Ylab-fParCollBaryonProj.GetPt2NuclearDestructP4()) );
753  coeff *= exfactor;
754  coeff /= ( 1. + exfactor );
755  SetPt2ofNuclearDestruction( (fParCollBaryonProj.GetPt2NuclearDestructP1()+coeff)*CLHEP::GeV*CLHEP::GeV );
756 
757  SetMaxPt2ofNuclearDestruction( fParCollBaryonProj.GetMaxPt2ofNuclearDestruct() );
758  SetExcitationEnergyPerWoundedNucleon( fParCollBaryonProj.GetExciEnergyPerWoundedNucleon() );
759 
760  }
761 
762  #ifdef debugFTFparams
763  G4cout<<"CofNuclearDestructionPr "<< GetCofNuclearDestructionPr() << G4endl;
764  G4cout<<"CofNuclearDestructionTr "<< GetCofNuclearDestruction() << G4endl;
765  G4cout<<"R2ofNuclearDestruction "<< GetR2ofNuclearDestruction()/fermi/fermi <<" fermi^2"<< G4endl;
766  G4cout<<"DofNuclearDestruction "<< GetDofNuclearDestruction() << G4endl;
767  G4cout<<"Pt2ofNuclearDestruction "<< GetPt2ofNuclearDestruction()/GeV/GeV <<" GeV^2"<< G4endl;
768  G4cout<<"ExcitationEnergyPerWoundedNucleon "<< GetExcitationEnergyPerWoundedNucleon() <<" MeV"<< G4endl;
769  #endif
770 
771  //SetCofNuclearDestruction( 0.47*G4Exp( 2.0*(Ylab - 2.5) )/( 1.0 + G4Exp( 2.0*(Ylab - 2.5) ) ) );
772  //SetPt2ofNuclearDestruction( ( 0.035 + 0.1*G4Exp( 4.0*(Ylab - 3.0) )/( 1.0 + G4Exp( 4.0*(Ylab - 3.0) ) ) )*GeV*GeV );
773 
774  //SetMagQuarkExchange( 120.0 ); // 210.0 PipP
775  //SetSlopeQuarkExchange( 2.0 );
776  //SetDeltaProbAtQuarkExchange( 0.6 );
777  //SetProjMinDiffMass( 0.7 ); // GeV 1.1
778  //SetProjMinNonDiffMass( 0.7 ); // GeV
779  //SetProbabilityOfProjDiff( 0.0); // 0.85*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
780  //SetTarMinDiffMass( 1.1 ); // GeV
781  //SetTarMinNonDiffMass( 1.1 ); // GeV
782  //SetProbabilityOfTarDiff( 0.0 ); // 0.85*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
783 
784  //SetAveragePt2( 0.0 ); // GeV^2 0.3
785  //------------------------------------
786  //SetProbabilityOfElasticScatt( 1.0, 1.0); //(Xtotal, Xelastic);
787  //SetProbabilityOfProjDiff( 1.0*0.62*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.51 ) ); // 0->1
788  //SetProbabilityOfTarDiff( 4.0*0.62*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.51 ) ); // 2->4
789  //SetAveragePt2( 0.3 ); // (0.15)
790  //SetAvaragePt2ofElasticScattering( 0.0 );
791 
792  //SetMaxNumberOfCollisions( Plab, 6.0 ); //(4.0*(Plab + 0.01), Plab); // 6.0 );
793  //SetAveragePt2( 0.15 );
794  //SetCofNuclearDestruction(-1.);//( 0.75 ); // (0.25)
795  //SetExcitationEnergyPerWoundedNucleon(0.);//( 30.0*MeV ); // (75.0*MeV)
796  //SetDofNuclearDestruction(0.);//( 0.2 ); //0.4 // 0.3 0.5
797 
798  /*
799  SetAveragePt2(0.3);
800  SetCofNuclearDestructionPr(0.);
801  SetCofNuclearDestruction(0.); //( 0.5 ); (0.25)
802  SetExcitationEnergyPerWoundedNucleon(0.); // 30.0*MeV; (75.0*MeV)
803  SetDofNuclearDestruction(0.); // 0.2; 0.4; 0.3; 0.5
804  SetPt2ofNuclearDestruction(0.); //(2.*0.075*GeV*GeV); ( 0.3*GeV*GeV ); (0.168*GeV*GeV)
805  */
806 
807  //SetExcitationEnergyPerWoundedNucleon(0.001);
808  //SetPt2Kink( 0.0*GeV*GeV );
809 
810  //SetRadiusOfHNinteractions2( Xtotal/pi/10.0 /2.);
811  //SetRadiusOfHNinteractions2( (Xtotal - Xelastic)/pi/10.0 );
812  //SetProbabilityOfElasticScatt( 1.0, 0.0);
813  /*
814  G4cout << "Pt2 " << GetAveragePt2()<<" "<<GetAveragePt2()/GeV/GeV<<G4endl;
815  G4cout << "Cnd " << GetCofNuclearDestruction() << G4endl;
816  G4cout << "Dnd " << GetDofNuclearDestruction() << G4endl;
817  G4cout << "Pt2 " << GetPt2ofNuclearDestruction()/GeV/GeV << G4endl;
818  */
819 
820 }
821 
822 //============================================================================
823 
824 G4double G4FTFParameters::GetMinMass( const G4ParticleDefinition* aParticle ) {
825  // The code is used for estimating the minimal string mass produced in diffraction dissociation.
826  // The indices used for minMassQDiQStr must be between 1 and 5, corresponding to the 5 considered
827  // quarks: d, u, s, c and b; enforcing this explicitly avoids compilation errors.
828  G4double EstimatedMass = 0.0;
829  G4int partID = std::abs(aParticle->GetPDGEncoding());
830  G4int Qleft = std::max( partID/100, 1 );
831  G4int Qright = std::max( (partID/ 10)%10, 1 );
832  if ( Qleft < 6 && Qright < 6 ) { // Q-Qbar string
833  EstimatedMass = StringMass->minMassQQbarStr[Qleft-1][Qright-1];
834  } else if ( Qleft < 6 && Qright > 6 ) { // Q - DiQ string
835  G4int q1 = std::max( std::min( Qright/10, 5 ), 1 );
836  G4int q2 = std::max( std::min( Qright%10, 5 ), 1 );
837  EstimatedMass = StringMass->minMassQDiQStr[Qleft-1][q1-1][q2-1];
838  } else if ( Qleft > 6 && Qright < 6 ) { // DiQ - Q string
839  G4int q1 = std::max( std::min( Qleft/10, 5 ), 1 );
840  G4int q2 = std::max( std::min( Qleft%10, 5 ), 1 );
841  EstimatedMass = StringMass->minMassQDiQStr[Qright-1][q1-1][q2-1];
842  }
843  return EstimatedMass;
844 }
845 
846 //============================================================================
847 
848 G4double G4FTFParameters::GetProcProb( const G4int ProcN, const G4double y ) {
849  G4double Prob( 0.0 );
850  if ( y < ProcParams[ProcN][6] ) {
851  Prob = ProcParams[ProcN][5];
852  if (Prob < 0.) Prob=0.;
853  return Prob;
854  }
855  Prob = ProcParams[ProcN][0] * G4Exp( -ProcParams[ProcN][1]*y ) +
856  ProcParams[ProcN][2] * G4Exp( -ProcParams[ProcN][3]*y ) +
857  ProcParams[ProcN][4];
858  if (Prob < 0.) Prob=0.;
859  return Prob;
860 }
861 
862 //============================================================================
863 
864 G4HadronicDeveloperParameters& HDP = G4HadronicDeveloperParameters::GetInstance();
865 
866 
867 class G4FTFSettingDefaultHDP
868 {
869  public:
870 
871  // ctor
872  G4FTFSettingDefaultHDP() {
873 
874  //==================================================================
875  //
876  // Cross sections for elementary processes
877  //
878  // these are for Inelastic interactions, i.e. Xinelastic=(Xtotal-Xelastix)>0.
879  // for elastic, all the A's & B's, Atop & Ymin are zeros
880  // general formula: Pp = A1*exp(B1*Y) + A2*exp(B2*Y) + A3
881  // but if Y<Ymin, then Pp=max(0.,Atop)
882  // for details, see also G4FTFParameters::GetProcProb( ProcN, y )
883  //
884  // baryons
885  //
886  /* JVY, Oct. 31, 2017: Per Alberto R. & Vladimir U., keep this group of parameters FIXED
887  // Process=0 --> Qexchg w/o excitation
888  //
889  HDP.SetDefault( "FTF_BARYON_PROC0_A1", 13.71 );
890  HDP.SetDefault( "FTF_BARYON_PROC0_B1", 1.75 );
891  HDP.SetDefault( "FTF_BARYON_PROC0_A2", -30.69 );
892  HDP.SetDefault( "FTF_BARYON_PROC0_B2", 3.0 );
893  HDP.SetDefault( "FTF_BARYON_PROC0_A3", 0.0 );
894  HDP.SetDefault( "FTF_BARYON_PROC0_ATOP", 1.0 );
895  HDP.SetDefault( "FTF_BARYON_PROC0_YMIN", 0.93 );
896  //
897  // Process=1 --> Qexchg w/excitation
898  //
899  HDP.SetDefault( "FTF_BARYON_PROC1_A1", 25. );
900  HDP.SetDefault( "FTF_BARYON_PROC1_B1", 1. );
901  HDP.SetDefault( "FTF_BARYON_PROC1_A2", -50.34 );
902  HDP.SetDefault( "FTF_BARYON_PROC1_B2", 1.5 );
903  HDP.SetDefault( "FTF_BARYON_PROC1_A3", 0. );
904  HDP.SetDefault( "FTF_BARYON_PROC1_ATOP", 0. );
905  HDP.SetDefault( "FTF_BARYON_PROC1_YMIN", 1.4 );
906  */
907  //
908  // NOTE: Process #2 & 3 are projectile & target diffraction
909  // they have more complex definition of A1 & A2
910  // (see example below)
911  // SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93);// Projectile diffraction
912  // SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93);// Target diffraction
913  //
914  // Also, for ( AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10 )
915  // projectile and/or target diffraction (dissociation) may be switched ON/OFF
916  //
917  HDP.SetDefault( "FTF_BARYON_DIFF_DISSO_PROJ", false );
918  HDP.SetDefault( "FTF_BARYON_DIFF_DISSO_TGT", false ); // as in hadr-string-diff-V10-03-07
919  //
920  /* JVY, Oct. 31, 2017: Per Alberto R. & Vladimir U., keep this group of parameters FIXED
921  // Process=4 --> Qexchg w/additional multiplier in excitation
922  //
923  HDP.SetDefault( "FTF_BARYON_PROC4_A1", 0.6 );
924  HDP.SetDefault( "FTF_BARYON_PROC4_B1", 0. );
925  HDP.SetDefault( "FTF_BARYON_PROC4_A2", -1.2 );
926  HDP.SetDefault( "FTF_BARYON_PROC4_B2", 0.5 );
927  HDP.SetDefault( "FTF_BARYON_PROC4_A3", 0. );
928  HDP.SetDefault( "FTF_BARYON_PROC4_ATOP",0. );
929  HDP.SetDefault( "FTF_BARYON_PROC4_YMIN",1.4 );
930  */
931  //
932  // Parameters of participating hadron (baryon) excitation
933  //
934  HDP.SetDefault( "FTF_BARYON_DELTA_PROB_QEXCHG", 0. );
935  HDP.SetDefault( "FTF_BARYON_PROB_SAME_QEXCHG", 0. );
936  HDP.SetDefault( "FTF_BARYON_DIFF_M_PROJ", 1.16, 1.16, 3. ); // it's supposed to be in GeV but do NOT do (*CLHEP::GeV)
937  // because it'll be done in the G4FTFParameters::SetProjMinDiffMass
938  HDP.SetDefault( "FTF_BARYON_NONDIFF_M_PROJ", 1.16, 1.16, 3. ); // do NOT (*CLHEP::GeV) - same as above
939  HDP.SetDefault( "FTF_BARYON_DIFF_M_TGT", 1.16, 1.16, 3. ); // do NOT (*CLHEP::GeV) - same as above
940  HDP.SetDefault( "FTF_BARYON_NONDIFF_M_TGT", 1.16, 1.16, 3. ); // do NOT (*CLHEP::GeV) - same as above
941  HDP.SetDefault( "FTF_BARYON_AVRG_PT2", 0.3, 0.08, 1. ); // do NOT (*CLHEP::GeV*CLHEP::GeV)
942  //
943  // JVY, Oct. 6, 2017: Per Alberto R., keep these two settings fixed (for now)
944  //
945  // HDP.SetDefault( "FTF_BARYON_PROB_DISTR_PROJ", 0.3 );
946  // HDP.SetDefault( "FTF_BARYON_PROB_DISTR_TGT", 0.3 );
947 
948 
949  // pions
950  //
951  // JVY, Aug.8, 2018 --> Feb.14, 2019 --> June 25, 2019:
952  // Parameters of participating hadron (pions) excitation
953  //
954  /* JVY, June 25, 2019: For now, keep this group of parameters FIXED
955  // Process=0 --> Qexchg w/o excitation
956  //
957  HDP.SetDefault( "FTF_PION_PROC0_A1", 150.0 );
958  HDP.SetDefault( "FTF_PION_PROC0_B1", 1.8 );
959  HDP.SetDefault( "FTF_PION_PROC0_A2",-247.3 );
960  HDP.SetDefault( "FTF_PION_PROC0_B2", 2.3 );
961  HDP.SetDefault( "FTF_PION_PROC0_A3", 0.0 );
962  HDP.SetDefault( "FTF_PION_PROC0_ATOP", 1.0 );
963  HDP.SetDefault( "FTF_PION_PROC0_YMIN", 2.3 );
964  //
965  // Process=1 --> Qexchg w/excitation
966  //
967  HDP.SetDefault( "FTF_PION_PROC1_A1", 5.77 );
968  HDP.SetDefault( "FTF_PION_PROC1_B1", 0.6 );
969  HDP.SetDefault( "FTF_PION_PROC1_A2", -5.77 );
970  HDP.SetDefault( "FTF_PION_PROC1_B2", 0.8 );
971  HDP.SetDefault( "FTF_PION_PROC1_A3", 0. );
972  HDP.SetDefault( "FTF_PION_PROC1_ATOP", 0. );
973  HDP.SetDefault( "FTF_PION_PROC1_YMIN", 0.0 );
974  //
975  // NOTE: Process #2 & 3 are projectile & target diffraction
976  //
977  // Process=2 --> Projectile diffraction
978  //
979  HDP.SetDefault( "FTF_PION_PROC2_A1", 2.27 );
980  HDP.SetDefault( "FTF_PION_PROC2_B1", 0.5 );
981  HDP.SetDefault( "FTF_PION_PROC2_A2", -98052.0);
982  HDP.SetDefault( "FTF_PION_PROC2_B2", 4.0 );
983  HDP.SetDefault( "FTF_PION_PROC2_A3", 0. );
984  HDP.SetDefault( "FTF_PION_PROC2_ATOP", 0. );
985  HDP.SetDefault( "FTF_PION_PROC2_YMIN", 3.0 );
986  //
987  // Process=3 --> Target diffraction
988  //
989  HDP.SetDefault( "FTF_PION_PROC3_A1", 7.0 );
990  HDP.SetDefault( "FTF_PION_PROC3_B1", 0.9 );
991  HDP.SetDefault( "FTF_PION_PROC3_A2", -85.28);
992  HDP.SetDefault( "FTF_PION_PROC3_B2", 1.9 );
993  HDP.SetDefault( "FTF_PION_PROC3_A3", 0.08);
994  HDP.SetDefault( "FTF_PION_PROC3_ATOP", 0. );
995  HDP.SetDefault( "FTF_PION_PROC3_YMIN", 2.2 );
996  */
997  //
998  // projectile and/or target diffraction (dissociation) may be switched ON/OFF
999  //
1000  // NOTE: Both projectile and target diffraction are turned OFF if
1001  // a) Number of Target Nucleons > 10 (NumberOfTargetNucleons > 10)
1002  // OR
1003  // b) Projectile Baryon Number > 10 (AbsProjectileBaryonNumber > 10)
1004  // ... which is "strange" because projectile is a pion !!!... so it's always OFF
1005  //
1006  HDP.SetDefault( "FTF_PION_DIFF_DISSO_PROJ", false );
1007  HDP.SetDefault( "FTF_PION_DIFF_DISSO_TGT", false );
1008  //
1009  /* JVY, June 25, 2019: For now keep this group of parameters FIXED
1010  // Process=4 --> Qexchg w/additional multiplier in excitation
1011  //
1012  HDP.SetDefault( "FTF_PION_PROC4_A1", 1.0 );
1013  HDP.SetDefault( "FTF_PION_PROC4_B1", 0. );
1014  HDP.SetDefault( "FTF_PION_PROC4_A2",-11.02);
1015  HDP.SetDefault( "FTF_PION_PROC4_B2", 1.0 );
1016  HDP.SetDefault( "FTF_PION_PROC4_A3", 0. );
1017  HDP.SetDefault( "FTF_PION_PROC4_ATOP",0. );
1018  HDP.SetDefault( "FTF_PION_PROC4_YMIN",2.4 );
1019  */
1020  //
1021  // NOTE; As of geant4-10-05, all these settings beloe are correct
1022  // (and are the same as they were in 10.4.ref06)
1023 
1024  HDP.SetDefault( "FTF_PION_DELTA_PROB_QEXCHG", 0.56 ); // in the past used to be 0.35
1025  HDP.SetDefault( "FTF_PION_DIFF_M_PROJ", 1.0, 0.5, 3. ); // in the past used to be 0.5...
1026  HDP.SetDefault( "FTF_PION_NONDIFF_M_PROJ", 1.0, 0.5, 3. ); // so why not set 0.5 as a low limit ?...
1027  // ... or perhaps even lower ?
1028  HDP.SetDefault( "FTF_PION_DIFF_M_TGT", 1.16, 1.16, 3. ); // All (NON)DIFF_M's are supposed to be in GeV but do NOT do (*CLHEP::GeV)
1029  // because it'll be done in the G4FTFParameters::SetProjMinDiffMass
1030  HDP.SetDefault( "FTF_PION_NONDIFF_M_TGT", 1.16, 1.16, 3. );
1031  HDP.SetDefault( "FTF_PION_AVRG_PT2", 0.3, 0.08, 1. ); // do NOT (*CLHEP::GeV*CLHEP::GeV)
1032 
1033  //==================================================================
1034  //
1035  // nuclear destruction
1036  //
1037  // NOTE: Settings of most of these parameters are the same
1038  // for different types of projectile hadron
1039  // However, we decided to introduce separate variables
1040  // and configuration cards for each type of projectile
1041  //
1042  // baryons
1043  //
1044  // projectile destruction
1045  //
1046  HDP.SetDefault( "FTF_BARYON_NUCDESTR_P1_PROJ", 1., 0., 1. ); // in principle, it should be 1./NBRN - FIXME later !
1047  HDP.SetDefault( "FTF_BARYON_NUCDESTR_P1_NBRN_PROJ", false );
1048  //
1049  // for now, keep fixed p2 & p3 for the proj destruction
1050  // they're defined explicitly in G4FTFParamCollection ctor
1051  //
1052  // target destruction
1053  //
1054  HDP.SetDefault( "FTF_BARYON_NUCDESTR_P1_TGT", 1., 0., 1. );
1055  HDP.SetDefault( "FTF_BARYON_NUCDESTR_P1_ADEP_TGT", false );
1056  HDP.SetDefault( "FTF_BARYON_NUCDESTR_P2_TGT", 4.0, 2., 16. );
1057  HDP.SetDefault( "FTF_BARYON_NUCDESTR_P3_TGT", 2.1, 0., 4. );
1058  //
1059  HDP.SetDefault( "FTF_BARYON_PT2_NUCDESTR_P1", 0.035, 0., 0.25 );
1060  HDP.SetDefault( "FTF_BARYON_PT2_NUCDESTR_P2", 0.04, 0., 0.25 );
1061  HDP.SetDefault( "FTF_BARYON_PT2_NUCDESTR_P3", 4.0, 2., 16. );
1062  HDP.SetDefault( "FTF_BARYON_PT2_NUCDESTR_P4", 2.5, 0., 4. );
1063  //
1064  HDP.SetDefault( "FTF_BARYON_NUCDESTR_R2", 1.5*CLHEP::fermi*CLHEP::fermi, 0.5*CLHEP::fermi*CLHEP::fermi, 2.*CLHEP::fermi*CLHEP::fermi );
1065  HDP.SetDefault( "FTF_BARYON_EXCI_E_PER_WNDNUCLN", 40.*CLHEP::MeV, 0., 100.*CLHEP::MeV );
1066  HDP.SetDefault( "FTF_BARYON_NUCDESTR_DISP", 0.3, 0.1, 0.4 );
1067  //
1068  // JVY, Oct. 6, 2017: Per Alberto R., this is just a technical parameter,
1069  // and it should NOT be changed
1070  //
1071  // HDP.SetDefault( "FTF_BARYON_NUCDESTR_MAXPT2", 1. * CLHEP::GeV*CLHEP::GeV );
1072 
1073 
1074  // mesons - these parameters are common for pions, kaons, etc. (per original code)
1075  //
1076  // NOTE: *NO* projectile destruction for mesons !!!
1077  //
1078  // target destruction
1079  //
1080  HDP.SetDefault( "FTF_MESON_NUCDESTR_P1_TGT", 0.00481, 0., 1. );
1081  HDP.SetDefault( "FTF_MESON_NUCDESTR_P1_ADEP_TGT", true );
1082  HDP.SetDefault( "FTF_MESON_NUCDESTR_P2_TGT", 4.0, 2., 16. );
1083  HDP.SetDefault( "FTF_MESON_NUCDESTR_P3_TGT", 2.1, 0., 4. );
1084  //
1085  HDP.SetDefault( "FTF_MESON_PT2_NUCDESTR_P1", 0.035, 0., 0.25 );
1086  HDP.SetDefault( "FTF_MESON_PT2_NUCDESTR_P2", 0.04, 0., 0.25 );
1087  HDP.SetDefault( "FTF_MESON_PT2_NUCDESTR_P3", 4.0, 2., 16. );
1088  HDP.SetDefault( "FTF_MESON_PT2_NUCDESTR_P4", 2.5, 0., 4. );
1089  //
1090  HDP.SetDefault( "FTF_MESON_NUCDESTR_R2", 1.5*CLHEP::fermi*CLHEP::fermi,
1091  0.5*CLHEP::fermi*CLHEP::fermi,
1092  2.*CLHEP::fermi*CLHEP::fermi );
1093  HDP.SetDefault( "FTF_MESON_EXCI_E_PER_WNDNUCLN", 40.*CLHEP::MeV, 0., 100.*CLHEP::MeV );
1094  HDP.SetDefault( "FTF_MESON_NUCDESTR_DISP", 0.3, 0.1, 0.4 );
1095 
1096  }
1097 };
1098 
1099 
1100 G4FTFSettingDefaultHDP FTFDefaultsHDP;
1101 
1102 
1103 G4FTFParamCollection::G4FTFParamCollection()
1104 {
1105 
1106  Reset(); // zero out everything
1107 
1108  //
1109  // keep the 2 parameters below fixed for now (i.e. do not take them from HDP)
1110  //
1111  fNuclearProjDestructP2 = 4.0;
1112  fNuclearProjDestructP3 = 2.1;
1113 
1114 }
1115 
1116 
1117 void G4FTFParamCollection::Reset()
1118 {
1119  // parameters of excitation
1120 
1121  // Proc=0 --> Qexchg w/o excitation
1122  fProc0A1 = 0.;
1123  fProc0B1 = 0.;
1124  fProc0A2 = 0.;
1125  fProc0B2 = 0.;
1126  fProc0A3 = 0.;
1127  fProc0Atop = 0.;
1128  fProc0Ymin = 0.;
1129 
1130  // Proc=1 --> Qexchg w/excitation
1131  fProc1A1 = 0.;
1132  fProc1B1 = 0.;
1133  fProc1A2 = 0.;
1134  fProc1B2 = 0.;
1135  fProc1A3 = 0.;
1136  fProc1Atop = 0.;
1137  fProc1Ymin = 0.;
1138 
1139  // Proc=2 & Proc=3 for ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 )
1140  // Do NOT do anything as it's set once and for all !!!
1141 
1142  fProjDiffDissociation = false;
1143  fTgtDiffDissociation = false;
1144 
1145  // Proc=4 --> Qexchg w/additional multiplier in excitation
1146  fProc4A1 = 0.;
1147  fProc4B1 = 0.;
1148  fProc4A2 = 0.;
1149  fProc4B2 = 0.;
1150  fProc4A3 = 0.;
1151  fProc4Atop = 0.;
1152  fProc4Ymin = 0.;
1153 
1154  // parameters of participating baryon excitation
1155 
1156  fDeltaProbAtQuarkExchange = 0.;
1157  fProbOfSameQuarkExchange = 0.;
1158  fProjMinDiffMass = 0.;
1159  fProjMinNonDiffMass = 0.;
1160  fTgtMinDiffMass = 0.;
1161  fTgtMinNonDiffMass = 0.;
1162  fAveragePt2 = 0.;
1163  fProbLogDistrPrD = 0.;
1164  fProbLogDistr = 0.;
1165 
1166  // parameters of nuclear distruction
1167 
1168  // COMMONs
1169  fNuclearProjDestructP1 = 0.;
1170  fNuclearTgtDestructP1 = 0.;
1171  fNuclearProjDestructP2 = 0.;
1172  fNuclearProjDestructP3 = 0.;
1173  fNuclearTgtDestructP2 = 0.;
1174  fNuclearTgtDestructP3 = 0.;
1175  fPt2NuclearDestructP1 = 0.;
1176  fPt2NuclearDestructP2 = 0.;
1177  fPt2NuclearDestructP3 = 0.;
1178  fPt2NuclearDestructP4 = 0.;
1179 
1180  // baryons
1181  fR2ofNuclearDestruct = 0.;
1182  fExciEnergyPerWoundedNucleon = 0.;
1183  fDofNuclearDestruct = 0.;
1184  fMaxPt2ofNuclearDestruct = 0.;
1185 
1186  return;
1187 }
1188 
1189 //============================================================================
1190 
1191 G4FTFParamCollBaryonProj::G4FTFParamCollBaryonProj()
1192  : G4FTFParamCollection()
1193 {
1194 
1195  // parameters of participating hadron (baryon) excitation
1196  //
1197  // baryons projectile
1198  //
1199  // Proc=0 --> Qexchg w/o excitation
1200  //
1201  /* As of Oct. 31, 2017 keep these fixed
1202  HDP.DeveloperGet( "FTF_BARYON_PROC0_A1", fProc0A1 );
1203  HDP.DeveloperGet( "FTF_BARYON_PROC0_B1", fProc0B1 );
1204  HDP.DeveloperGet( "FTF_BARYON_PROC0_A2", fProc0A2 );
1205  HDP.DeveloperGet( "FTF_BARYON_PROC0_B2", fProc0B2 );
1206  HDP.DeveloperGet( "FTF_BARYON_PROC0_A3", fProc0A3 );
1207  HDP.DeveloperGet( "FTF_BARYON_PROC0_ATOP", fProc0Atop );
1208  HDP.DeveloperGet( "FTF_BARYON_PROC0_YMIN", fProc0Ymin );
1209  */
1210  //
1211  fProc0A1 = 13.71;
1212  fProc0B1 = 1.75;
1213  fProc0A2 = -30.69; // (or -214.5 as in Doc ?)
1214  fProc0B2 = 3.; // ( or 4. as in Doc ?)
1215  fProc0A3 = 0.;
1216  fProc0Atop = 1.; // ( or 0.5 as in Doc ?)
1217  fProc0Ymin = 0.93; // (or 1.1 as in Doc ?)
1218  //
1219  // Proc=1 --> Qexchg w/excitation
1220  //
1221  /* As of Oct. 31, 2017 keep these fixed
1222  HDP.DeveloperGet( "FTF_BARYON_PROC1_A1", fProc1A1 );
1223  HDP.DeveloperGet( "FTF_BARYON_PROC1_B1", fProc1B1 );
1224  HDP.DeveloperGet( "FTF_BARYON_PROC1_A2", fProc1A2 );
1225  HDP.DeveloperGet( "FTF_BARYON_PROC1_B2", fProc1B2 );
1226  HDP.DeveloperGet( "FTF_BARYON_PROC1_A3", fProc1A3 );
1227  HDP.DeveloperGet( "FTF_BARYON_PROC1_ATOP", fProc1Atop );
1228  HDP.DeveloperGet( "FTF_BARYON_PROC1_YMIN", fProc1Ymin );
1229  */
1230  //
1231  fProc1A1 = 25.;
1232  fProc1B1 = 1.;
1233  fProc1A2 = -50.34;
1234  fProc1B2 = 1.5;
1235  fProc1A3 = 0.;
1236  fProc1Atop = 0.;
1237  fProc1Ymin = 1.4;
1238  //
1239  // Proc=2 & Proc=3 for the case ( AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10 )
1240  // (diffraction dissociation)
1241  // NOTE-1: used to be ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 )...
1242  // NOTE-2: As of 10.5, both are set to false (via HDP)
1243  //
1244  HDP.DeveloperGet( "FTF_BARYON_DIFF_DISSO_PROJ", fProjDiffDissociation );
1245  HDP.DeveloperGet( "FTF_BARYON_DIFF_DISSO_TGT", fTgtDiffDissociation );
1246  //
1247  //
1248  // Proc=4 --> Qexchg "w/additional multiplier" in excitation
1249  //
1250  /* As of Oct. 31, 2017 keep these fixed
1251  HDP.DeveloperGet( "FTF_BARYON_PROC4_A1", fProc4A1 );
1252  HDP.DeveloperGet( "FTF_BARYON_PROC4_B1", fProc4B1 );
1253  HDP.DeveloperGet( "FTF_BARYON_PROC4_A2", fProc4A2 );
1254  HDP.DeveloperGet( "FTF_BARYON_PROC4_B2", fProc4B2 );
1255  HDP.DeveloperGet( "FTF_BARYON_PROC4_A3", fProc4A3 );
1256  HDP.DeveloperGet( "FTF_BARYON_PROC4_ATOP", fProc4Atop );
1257  HDP.DeveloperGet( "FTF_BARYON_PROC4_YMIN", fProc4Ymin );
1258  */
1259  //
1260  fProc4A1 = 0.6; // (or 1. as in Doc ?)
1261  fProc4B1 = 0.;
1262  fProc4A2 = -1.2; // (or -2.01 as in Doc ?)
1263  fProc4B2 = 0.5;
1264  fProc4A3 = 0.;
1265  fProc4Atop = 0.;
1266  fProc4Ymin = 1.4;
1267  //
1268  //
1269  HDP.DeveloperGet( "FTF_BARYON_DELTA_PROB_QEXCHG", fDeltaProbAtQuarkExchange );
1270  HDP.DeveloperGet( "FTF_BARYON_PROB_SAME_QEXCHG", fProbOfSameQuarkExchange );
1271  HDP.DeveloperGet( "FTF_BARYON_DIFF_M_PROJ", fProjMinDiffMass );
1272  HDP.DeveloperGet( "FTF_BARYON_NONDIFF_M_PROJ", fProjMinNonDiffMass );
1273  HDP.DeveloperGet( "FTF_BARYON_DIFF_M_TGT", fTgtMinDiffMass );
1274  HDP.DeveloperGet( "FTF_BARYON_NONDIFF_M_TGT", fTgtMinNonDiffMass );
1275  HDP.DeveloperGet( "FTF_BARYON_AVRG_PT2", fAveragePt2 );
1276  //
1277  // JVY - Per Alberto R., we're curretly keeping these two settings fixed,
1278  // thus they're defined here explicitly, rather than via HDP
1279  //
1280  // HDP.DeveloperGet( "FTF_BARYON_PROB_DISTR_PROJ", fProbLogDistrPrD );
1281  // HDP.DeveloperGet( "FTF_BARYON_PROB_DISTR_TGT", fProbLogDistr );
1282  fProbLogDistrPrD = 0.55;
1283  fProbLogDistr = 0.55;
1284 
1285  // nuclear destruction
1286  //
1287  // ---> LATER !!! ---> fBaryonMaxNumberOfCollisions = 2.;
1288  //
1289 
1290  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_P1_PROJ", fNuclearProjDestructP1 );
1291  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_P1_NBRN_PROJ",fNuclearProjDestructP1_NBRNDEP );
1292  //
1293  // keep the 2 parameters below fixed for now (i.e. do not take them from HDP)
1294  //
1295  fNuclearProjDestructP2 = 4.0;
1296  fNuclearProjDestructP3 = 2.1;
1297 
1298  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_P1_TGT", fNuclearTgtDestructP1 );
1299  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_P1_ADEP_TGT", fNuclearTgtDestructP1_ADEP );
1300  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_P2_TGT", fNuclearTgtDestructP2 );
1301  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_P3_TGT", fNuclearTgtDestructP3 );
1302  //
1303  HDP.DeveloperGet( "FTF_BARYON_PT2_NUCDESTR_P1", fPt2NuclearDestructP1 );
1304  HDP.DeveloperGet( "FTF_BARYON_PT2_NUCDESTR_P2", fPt2NuclearDestructP2 );
1305  HDP.DeveloperGet( "FTF_BARYON_PT2_NUCDESTR_P3", fPt2NuclearDestructP3 );
1306  HDP.DeveloperGet( "FTF_BARYON_PT2_NUCDESTR_P4", fPt2NuclearDestructP4 );
1307  //
1308  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_R2", fR2ofNuclearDestruct );
1309  HDP.DeveloperGet( "FTF_BARYON_EXCI_E_PER_WNDNUCLN", fExciEnergyPerWoundedNucleon );
1310  //
1311  HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_DISP", fDofNuclearDestruct ); // NOTE: "Dof" means "Dispersion of..."
1312  //
1313  // NOTE-1: this parameter has changed from 1. to 9. between 10.2 and 10.3.ref07 !!!
1314  // ... then it went back to 1. for the 10.4-candidate...
1315  // NOTE-2: this is a "technical" parameter, it should not be changed; this is why
1316  // it is defined explicitly rather than via HDP
1317  // --> HDP.DeveloperGet( "FTF_BARYON_NUCDESTR_MAXPT2", fMaxPt2ofNuclearDestruct );
1318  fMaxPt2ofNuclearDestruct = 9. * CLHEP::GeV*CLHEP::GeV;
1319 }
1320 
1321 
1322 G4FTFParamCollMesonProj::G4FTFParamCollMesonProj()
1323  : G4FTFParamCollection()
1324 {
1325  // nuclear destruction
1326  //
1327  // JVY, Aug.8, 2018 --> Feb.14, 2018 --> June 25, 2019:
1328  // These parameters are common for all mesons
1329  //
1330 
1331  HDP.DeveloperGet( "FTF_MESON_NUCDESTR_P1_TGT", fNuclearTgtDestructP1 );
1332  HDP.DeveloperGet( "FTF_MESON_NUCDESTR_P1_ADEP_TGT", fNuclearTgtDestructP1_ADEP );
1333  HDP.DeveloperGet( "FTF_MESON_NUCDESTR_P2_TGT", fNuclearTgtDestructP2 );
1334  HDP.DeveloperGet( "FTF_MESON_NUCDESTR_P3_TGT", fNuclearTgtDestructP3 );
1335  //
1336  HDP.DeveloperGet( "FTF_MESON_PT2_NUCDESTR_P1", fPt2NuclearDestructP1 );
1337  HDP.DeveloperGet( "FTF_MESON_PT2_NUCDESTR_P2", fPt2NuclearDestructP2 );
1338  HDP.DeveloperGet( "FTF_MESON_PT2_NUCDESTR_P3", fPt2NuclearDestructP3 );
1339  HDP.DeveloperGet( "FTF_MESON_PT2_NUCDESTR_P4", fPt2NuclearDestructP4 );
1340  //
1341  HDP.DeveloperGet( "FTF_MESON_NUCDESTR_R2", fR2ofNuclearDestruct );
1342  HDP.DeveloperGet( "FTF_MESON_EXCI_E_PER_WNDNUCLN", fExciEnergyPerWoundedNucleon );
1343  HDP.DeveloperGet( "FTF_MESON_NUCDESTR_DISP", fDofNuclearDestruct ); // NOTE: "Dof" means "Dispersion of..."
1344  //
1345  // NOTE: it is a "technical" parameter, it should not be changed;
1346  // this is why it is defined explicitly rather than via HDP
1347  //
1348  fMaxPt2ofNuclearDestruct = 1. * CLHEP::GeV*CLHEP::GeV;
1349 }
1350 
1351 
1352 G4FTFParamCollPionProj::G4FTFParamCollPionProj()
1353  : G4FTFParamCollMesonProj()
1354 {
1355  // parameters of participating pion excitation (pi+/- or pi0)
1356  //
1357  // Proc=0 --> Qexchg w/o excitation
1358  //
1359  /* As of June 25, 2019 keep these fixed
1360  HDP.DeveloperGet( "FTF_PION_PROC0_A1", fProc0A1 );
1361  HDP.DeveloperGet( "FTF_PION_PROC0_B1", fProc0B1 );
1362  HDP.DeveloperGet( "FTF_PION_PROC0_A2", fProc0A2 );
1363  HDP.DeveloperGet( "FTF_PION_PROC0_B2", fProc0B2 );
1364  HDP.DeveloperGet( "FTF_PION_PROC0_A3", fProc0A3 );
1365  HDP.DeveloperGet( "FTF_PION_PROC0_ATOP", fProc0Atop );
1366  HDP.DeveloperGet( "FTF_PION_PROC0_YMIN", fProc0Ymin );
1367  */
1368  //
1369  fProc0A1 = 150.0;
1370  fProc0B1 = 1.8;
1371  fProc0A2 =-247.3;
1372  fProc0B2 = 2.3; // ( or 4. as in Doc ?)
1373  fProc0A3 = 0.;
1374  fProc0Atop = 1.; // ( or 0.5 as in Doc ?)
1375  fProc0Ymin = 2.3; // (or 1.1 as in Doc ?)
1376  //
1377  // Proc=1 --> Qexchg w/excitation
1378  //
1379  /* As of Oct. 31, 2017 keep these fixed
1380  HDP.DeveloperGet( "FTF_PION_PROC1_A1", fProc1A1 );
1381  HDP.DeveloperGet( "FTF_PION_PROC1_B1", fProc1B1 );
1382  HDP.DeveloperGet( "FTF_PION_PROC1_A2", fProc1A2 );
1383  HDP.DeveloperGet( "FTF_PION_PROC1_B2", fProc1B2 );
1384  HDP.DeveloperGet( "FTF_PION_PROC1_A3", fProc1A3 );
1385  HDP.DeveloperGet( "FTF_PION_PROC1_ATOP", fProc1Atop );
1386  HDP.DeveloperGet( "FTF_PION_PROC1_YMIN", fProc1Ymin );
1387  */
1388  //
1389  fProc1A1 = 5.77;
1390  fProc1B1 = 0.6;
1391  fProc1A2 = -5.77;
1392  fProc1B2 = 0.8;
1393  fProc1A3 = 0.;
1394  fProc1Atop = 0.;
1395  fProc1Ymin = 0.0;
1396  //
1397  // Proc=2 --> Projectile diffraction
1398  //
1399  /* As of Oct. 31, 2017 keep these fixed
1400  HDP.DeveloperGet( "FTF_PION_PROC2_A1", fProc2A1 );
1401  HDP.DeveloperGet( "FTF_PION_PROC2_B1", fProc2B1 );
1402  HDP.DeveloperGet( "FTF_PION_PROC2_A2", fProc2A2 );
1403  HDP.DeveloperGet( "FTF_PION_PROC2_B2", fProc2B2 );
1404  HDP.DeveloperGet( "FTF_PION_PROC2_A3", fProc2A3 );
1405  HDP.DeveloperGet( "FTF_PION_PROC2_ATOP", fProc2Atop );
1406  HDP.DeveloperGet( "FTF_PION_PROC2_YMIN", fProc2Ymin );
1407  */
1408  //
1409  fProc2A1 = 2.27;
1410  fProc2B1 = 0.5;
1411  fProc2A2 =-98052.0;
1412  fProc2B2 = 4.0;
1413  fProc2A3 = 0.;
1414  fProc2Atop = 0.;
1415  fProc2Ymin = 3.0;
1416  //
1417  // Proc=3 --> Target diffraction
1418  //
1419  /* As of Oct. 31, 2017 keep these fixed
1420  HDP.DeveloperGet( "FTF_PION_PROC3_A1", fProc3A1 );
1421  HDP.DeveloperGet( "FTF_PION_PROC3_B1", fProc3B1 );
1422  HDP.DeveloperGet( "FTF_PION_PROC3_A2", fProc3A2 );
1423  HDP.DeveloperGet( "FTF_PION_PROC3_B2", fProc3B2 );
1424  HDP.DeveloperGet( "FTF_PION_PROC3_A3", fProc3A3 );
1425  HDP.DeveloperGet( "FTF_PION_PROC3_ATOP", fProc3Atop );
1426  HDP.DeveloperGet( "FTF_PION_PROC3_YMIN", fProc3Ymin );
1427  */
1428  //
1429  fProc3A1 = 7.0;
1430  fProc3B1 = 0.9;
1431  fProc3A2 = -85.28;
1432  fProc3B2 = 1.9;
1433  fProc3A3 = 0.08;
1434  fProc3Atop = 0.;
1435  fProc3Ymin = 2.2;
1436  //
1437  // for Proc2 & Proc3, pprojectile or target diffraction can be turned ON/OFF
1438  // if num.baryons >10 (which is strange for projectile which is pion !!!)
1439  //
1440  HDP.DeveloperGet( "FTF_PION_DIFF_DISSO_PROJ", fProjDiffDissociation );
1441  HDP.DeveloperGet( "FTF_PION_DIFF_DISSO_TGT", fTgtDiffDissociation );
1442  //
1443  // Proc=4 --> Qexchg "w/additional multiplier" in excitation
1444  //
1445  /* As of Oct. 31, 2017 keep these fixed
1446  HDP.DeveloperGet( "FTF_PION_PROC4_A1", fProc4A1 );
1447  HDP.DeveloperGet( "FTF_PION_PROC4_B1", fProc4B1 );
1448  HDP.DeveloperGet( "FTF_PION_PROC4_A2", fProc4A2 );
1449  HDP.DeveloperGet( "FTF_PION_PROC4_B2", fProc4B2 );
1450  HDP.DeveloperGet( "FTF_PION_PROC4_A3", fProc4A3 );
1451  HDP.DeveloperGet( "FTF_PION_PROC4_ATOP", fProc4Atop );
1452  HDP.DeveloperGet( "FTF_PION_PROC4_YMIN", fProc4Ymin );
1453  */
1454  //
1455  fProc4A1 = 1.0;
1456  fProc4B1 = 0.;
1457  fProc4A2 = -11.02;
1458  fProc4B2 = 1.0;
1459  fProc4A3 = 0.;
1460  fProc4Atop = 0.;
1461  fProc4Ymin = 2.4;
1462  //
1463  //
1464  HDP.DeveloperGet( "FTF_PION_DELTA_PROB_QEXCHG", fDeltaProbAtQuarkExchange );
1465  HDP.DeveloperGet( "FTF_PION_DIFF_M_PROJ", fProjMinDiffMass );
1466  HDP.DeveloperGet( "FTF_PION_NONDIFF_M_PROJ", fProjMinNonDiffMass );
1467  HDP.DeveloperGet( "FTF_PION_DIFF_M_TGT", fTgtMinDiffMass );
1468  HDP.DeveloperGet( "FTF_PION_NONDIFF_M_TGT", fTgtMinNonDiffMass );
1469  HDP.DeveloperGet( "FTF_PION_AVRG_PT2", fAveragePt2 );
1470  //
1471  fProbOfSameQuarkExchange = 0.; // This does NOT seem to apply to the pion case
1472  //
1473  // currently keep these two parameters fixed
1474  // thus they're defined here explicitly, rather than via HDP
1475  //
1476  fProbLogDistrPrD = 0.55;
1477  fProbLogDistr = 0.55;
1478 }
1479 
1480 
1481 //============================================================================
1482 
1483 G4FTFParameters::~G4FTFParameters() {
1484  if ( StringMass ) delete StringMass;
1485 }
1486 
1487 //============================================================================
1488 
1489 void G4FTFParameters::Reset()
1490 {
1491  FTFhNcmsEnergy = 0.0;
1492  FTFXtotal = 0.0;
1493  FTFXelastic = 0.0;
1494  FTFXinelastic = 0.0;
1495  FTFXannihilation = 0.0;
1496  ProbabilityOfAnnihilation = 0.0;
1497  ProbabilityOfElasticScatt = 0.0;
1498  RadiusOfHNinteractions2 = 0.0;
1499  FTFSlope = 0.0;
1500  AvaragePt2ofElasticScattering = 0.0;
1501  FTFGamma0 = 0.0;
1502  DeltaProbAtQuarkExchange = 0.0;
1503  ProbOfSameQuarkExchange = 0.0;
1504  ProjMinDiffMass = 0.0;
1505  ProjMinNonDiffMass = 0.0;
1506  ProbLogDistrPrD = 0.0;
1507  TarMinDiffMass = 0.0;
1508  TarMinNonDiffMass = 0.0;
1509  AveragePt2 = 0.0;
1510  ProbLogDistr = 0.0;
1511  Pt2kink = 0.0;
1512  MaxNumberOfCollisions = 0.0;
1513  ProbOfInelInteraction = 0.0;
1514  CofNuclearDestructionPr = 0.0;
1515  CofNuclearDestruction = 0.0;
1516  R2ofNuclearDestruction = 0.0;
1517  ExcitationEnergyPerWoundedNucleon = 0.0;
1518  DofNuclearDestruction = 0.0;
1519  Pt2ofNuclearDestruction = 0.0;
1520  MaxPt2ofNuclearDestruction = 0.0;
1521 
1522  for ( G4int i = 0; i < 4; i++ ) {
1523  for ( G4int j = 0; j < 7; j++ ) {
1524  ProcParams[i][j] = 0.0;
1525  }
1526  }
1527 
1528  return;
1529 }
1530