Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/diffraction/src/G4FTFParameters.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 #include "G4HadronicParameters.hh"
 56 
 57 //============================================================================
 58 
 59 //#define debugFTFparams
 60 
 61 //============================================================================
 62 
 63 G4FTFParameters::G4FTFParameters() 
 64 {
 65   // Set-up alternative sets of FTF parameters (called "tunes").
 66   // Note that the very first tune (with indexTune == 0) corresponds to the default
 67   // set of parameters, which does not need to be set-up explicitly: that's why
 68   // the for loop below starts from 1 and not from 0.
 69   // The check whether an alternative tune has been switched on is done at the
 70   // level of the G4FTFParamCollection::SetTune method.
 71   for ( G4int indexTune = 1; indexTune < G4FTFTunings::sNumberOfTunes; ++indexTune ) {
 72     fArrayParCollBaryonProj[indexTune].SetTune(indexTune);
 73     fArrayParCollMesonProj[indexTune].SetTune(indexTune);
 74     fArrayParCollPionProj[indexTune].SetTune(indexTune);
 75   }
 76   
 77   StringMass = new G4LundStringFragmentation;  // for estimation of min. mass of diffr. states
 78   Reset();
 79   csGGinstance = 
 80     G4CrossSectionDataSetRegistry::Instance()->GetComponentCrossSection("Glauber-Gribov");
 81   if (!csGGinstance) {
 82     csGGinstance = new G4ComponentGGHadronNucleusXsc();
 83   }
 84 
 85   EnableDiffDissociationForBGreater10 = G4HadronicParameters::Instance()->EnableDiffDissociationForBGreater10();
 86   
 87   // Set parameters of a string kink
 88   SetPt2Kink( 0.0*GeV*GeV );  // To switch off kinky strings (bad results obtained with 6.0*GeV*GeV)
 89   G4double Puubar( 1.0/3.0 ), Pddbar( 1.0/3.0 ), Pssbar( 1.0/3.0 );  // SU(3) symmetry
 90   //G4double Puubar( 0.41 ), Pddbar( 0.41 ), Pssbar( 0.18 );         // Broken SU(3) symmetry
 91   SetQuarkProbabilitiesAtGluonSplitUp( Puubar, Pddbar, Pssbar );
 92 }
 93 
 94 //============================================================================
 95 
 96 void G4FTFParameters::InitForInteraction( const G4ParticleDefinition* particle, 
 97                                           G4int theA, G4int theZ, G4double PlabPerParticle ) 
 98 {
 99   Reset();
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       = std::abs( 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( 0.0 ), Xelastic( 0.0 ), Xannihilation( 0.0 );
139   G4int NumberOfTargetNucleons;
140 
141   Ylab = 0.5 * G4Log( (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   Plab /= GeV;
154   G4double Xftf = 0.0;
155 
156   G4int NumberOfTargetProtons  = theZ; 
157   G4int NumberOfTargetNeutrons = theA - theZ;
158   NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
159 
160   // ---------- hadron projectile ----------------
161   if ( AbsProjectileBaryonNumber <= 1 ) {  // Projectile is hadron or baryon 
162 
163     // Interaction on P
164     G4double xTtP = csGGinstance->GetTotalIsotopeCrossSection( particle, KineticEnergy, 1, 1);
165     G4double xElP = csGGinstance->GetElasticIsotopeCrossSection(particle, KineticEnergy, 1, 1);
166 
167     // Interaction on N
168     G4double xTtN = csGGinstance->GetTotalIsotopeCrossSection( particle, KineticEnergy, 0, 1);
169     G4double xElN = csGGinstance->GetElasticIsotopeCrossSection(particle, KineticEnergy, 0, 1);
170 
171     // Average properties of h+N interactions
172     Xtotal   = ( NumberOfTargetProtons * xTtP + NumberOfTargetNeutrons * xTtN ) / NumberOfTargetNucleons;
173     Xelastic = ( NumberOfTargetProtons * xElP + NumberOfTargetNeutrons * xElN ) / NumberOfTargetNucleons; 
174     Xannihilation = 0.0;
175 
176     Xtotal /= millibarn;
177     Xelastic /= millibarn;
178 
179     #ifdef debugFTFparams
180     G4cout<<"Estimated cross sections (total and elastic) of h+N interactions "<<Xtotal<<" "<<Xelastic<<" (mb)"<<G4endl;
181     #endif
182   }
183 
184   // ---------- nucleus projectile ----------------
185   if ( ProjectileIsNucleus  &&  ProjectileBaryonNumber > 1 ) {
186 
187     #ifdef debugFTFparams
188     G4cout<<"Projectile is a nucleus: A and Z - "<<ProjectileBaryonNumber<<" "<<ProjectileCharge<<G4endl;
189     #endif
190 
191     const G4ParticleDefinition* Proton = G4Proton::Proton(); 
192     // Interaction on P
193     G4double XtotPP = csGGinstance->GetTotalIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
194     G4double XelPP  = csGGinstance->GetElasticIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
195 
196     const G4ParticleDefinition* Neutron = G4Neutron::Neutron();
197     // Interaction on N
198     G4double XtotPN = csGGinstance->GetTotalIsotopeCrossSection(Neutron, KineticEnergy, 0, 1);
199     G4double XelPN  = csGGinstance->GetElasticIsotopeCrossSection(Neutron, KineticEnergy, 0, 1);
200 
201     #ifdef debugFTFparams
202     G4cout << "XsPP (total and elastic) " << XtotPP/millibarn << " " << XelPP/millibarn <<" (mb)"<< G4endl
203            << "XsPN (total and elastic) " << XtotPN/millibarn << " " << XelPN/millibarn <<" (mb)"<< G4endl;
204     #endif
205 
206     Xtotal = ( 
207                   AbsProjectileCharge * NumberOfTargetProtons * XtotPP + 
208                   ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
209                       NumberOfTargetNeutrons * XtotPP 
210                 + 
211                   ( AbsProjectileCharge * NumberOfTargetNeutrons +
212                     ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
213                         NumberOfTargetProtons ) * XtotPN
214                 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
215     Xelastic= (
216                   AbsProjectileCharge * NumberOfTargetProtons * XelPP + 
217                   ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
218                       NumberOfTargetNeutrons * XelPP 
219                  + 
220                   ( AbsProjectileCharge * NumberOfTargetNeutrons +
221                     ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
222                         NumberOfTargetProtons ) * XelPN
223                 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
224 
225     Xannihilation = 0.0;
226     Xtotal /= millibarn;
227     Xelastic /= millibarn;
228   }
229 
230   // ---------- The projectile is anti-baryon or anti-nucleus ----------------
231   //                     anti  Sigma^0_c                  anti Delta^-
232   if ( ProjectilePDGcode >= -4112  &&  ProjectilePDGcode <= -1114 ) {
233     // Only non-strange and strange baryons are considered
234 
235     #ifdef debugFTFparams
236     G4cout<<"Projectile is a anti-baryon or anti-nucleus  - "<<ProjectileBaryonNumber<<" "<<ProjectileCharge<<G4endl;
237     G4cout<<"(Only non-strange and strange baryons are considered)"<<G4endl;
238     #endif
239 
240     G4double X_a( 0.0 ), X_b( 0.0 ), X_c( 0.0 ), X_d( 0.0 );
241     G4double MesonProdThreshold = ProjectileMass + TargetMass + 
242                                   ( 2.0 * 0.14 + 0.016 ); // 2 Mpi + DeltaE;
243 
244     if ( PlabPerParticle < 40.0*MeV ) { // Low energy limits. Projectile at rest.
245       Xtotal =   1512.9;    // mb
246       Xelastic =  473.2;    // mb
247       X_a =       625.1;    // mb
248       X_b =         9.780;  // mb
249       X_c =        49.989;  // mb
250       X_d =         6.614;  // mb
251     } else { // Total and elastic cross section of PbarP interactions a'la Arkhipov
252       G4double LogS = G4Log( ECMSsqr / 33.0625 );
253       G4double Xasmpt = 36.04 + 0.304*LogS*LogS;  // mb
254       LogS = G4Log( SqrtS / 20.74 );
255       G4double Basmpt = 11.92 + 0.3036*LogS*LogS;  // GeV^(-2)
256       G4double R0 = std::sqrt( 0.40874044*Xasmpt - Basmpt );  // GeV^(-1)
257 
258       G4double FlowF = SqrtS / std::sqrt( ECMSsqr*ECMSsqr + ProjectileMass2*ProjectileMass2 +
259                                           TargetMass2*TargetMass2 - 2.0*ECMSsqr*ProjectileMass2
260                                           - 2.0*ECMSsqr*TargetMass2 
261                                           - 2.0*ProjectileMass2*TargetMass2 );
262 
263       Xtotal = Xasmpt * ( 1.0 + 13.55*FlowF/R0/R0/R0*
264                                 (1.0 - 4.47/SqrtS + 12.38/ECMSsqr - 12.43/SqrtS/ECMSsqr) );  // mb
265 
266       Xasmpt = 4.4 + 0.101*LogS*LogS;  // mb
267       Xelastic = Xasmpt * ( 1.0 + 59.27*FlowF/R0/R0/R0*
268                                   (1.0 - 6.95/SqrtS + 23.54/ECMSsqr - 25.34/SqrtS/ECMSsqr ) ); // mb
269 
270       X_a = 25.0*FlowF;  // mb, 3-shirts diagram
271 
272       if ( SqrtS < MesonProdThreshold ) {
273         X_b = 3.13 + 140.0*G4Pow::GetInstance()->powA( MesonProdThreshold - SqrtS, 2.5 );  // mb anti-quark-quark annihilation
274         Xelastic -= 3.0*X_b;  // Xel-X(PbarP->NNbar)
275       } else {
276         X_b = 6.8/SqrtS;  // mb anti-quark-quark annihilation
277         Xelastic -= 3.0*X_b;  // Xel-X(PbarP->NNbar)
278       }
279 
280       X_c = 2.0*FlowF*sqr( ProjectileMass + TargetMass )/ECMSsqr;  // mb rearrangement
281 
282       X_d = 23.3/ECMSsqr;  // mb anti-quark-quark string creation
283     }
284 
285     G4double Xann_on_P( 0.0), Xann_on_N( 0.0 );
286 
287     if ( ProjectilePDGcode == -2212 ) {              // Pbar+P/N
288       Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0; 
289       Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
290     } else if ( ProjectilePDGcode == -2112 ) {       // NeutrBar+P/N
291       Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
292       Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0;
293     } else if ( ProjectilePDGcode == -3122 ) {       // LambdaBar+P/N
294       Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
295       Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
296     } else if ( ProjectilePDGcode == -3112 ) {       // Sigma-Bar+P/N
297       Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
298       Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
299     } else if ( ProjectilePDGcode == -3212 ) {       // Sigma0Bar+P/N
300       Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
301       Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
302     } else if ( ProjectilePDGcode == -3222 ) {       // Sigma+Bar+P/N
303       Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
304       Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
305     } else if ( ProjectilePDGcode == -3312 ) {       // Xi-Bar+P/N
306       Xann_on_P = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
307       Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
308     } else if ( ProjectilePDGcode == -3322 ) {       // Xi0Bar+P/N
309       Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
310       Xann_on_N = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
311     } else if ( ProjectilePDGcode == -3334 ) {       // Omega-Bar+P/N
312       Xann_on_P = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
313       Xann_on_N = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
314     } else {
315       G4cout << "Unknown anti-baryon for FTF annihilation" << G4endl;
316     }
317 
318     //G4cout << "Sum          " << Xann_on_P << G4endl;
319 
320     if ( ! ProjectileIsNucleus ) {  // Projectile is anti-baryon
321       Xannihilation = ( NumberOfTargetProtons * Xann_on_P  + NumberOfTargetNeutrons * Xann_on_N  )
322                       / NumberOfTargetNucleons;
323     } else {  // Projectile is a nucleus
324       Xannihilation = (
325                         ( AbsProjectileCharge * NumberOfTargetProtons + 
326                           ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
327                           NumberOfTargetNeutrons ) * Xann_on_P 
328                        + 
329                         ( AbsProjectileCharge * NumberOfTargetNeutrons +
330                           ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
331                           NumberOfTargetProtons ) * Xann_on_N
332                       ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
333     }
334 
335     //G4double Xftf = 0.0;  
336     MesonProdThreshold = ProjectileMass + TargetMass + (0.14 + 0.08); // Mpi + DeltaE
337     if ( SqrtS > MesonProdThreshold ) {
338       Xftf = 36.0 * ( 1.0 - MesonProdThreshold/SqrtS );
339     }
340 
341     Xtotal = Xelastic + Xannihilation + Xftf;
342 
343     #ifdef debugFTFparams
344     G4cout << "Plab Xtotal, Xelastic  Xinel Xftf " << Plab << " " << Xtotal << " " << Xelastic
345            << " " << Xtotal - Xelastic << " " << Xtotal - Xelastic - Xannihilation << " (mb)"<< G4endl
346            << "Plab Xelastic/Xtotal,  Xann/Xin " << Plab << " " << Xelastic/Xtotal << " " 
347            << Xannihilation/(Xtotal - Xelastic) << G4endl;
348     #endif
349 
350   }
351 
352   if ( Xtotal == 0.0 ) {  // Projectile is undefined, Nucleon assumed
353 
354     const G4ParticleDefinition* Proton = G4Proton::Proton(); 
355     // Interaction on P
356     G4double XtotPP = csGGinstance->GetTotalIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
357     G4double XelPP  = csGGinstance->GetElasticIsotopeCrossSection(Proton, KineticEnergy, 1, 1);
358 
359     // Interaction on N
360     G4double XtotPN = csGGinstance->GetTotalIsotopeCrossSection(Proton, KineticEnergy, 0, 1);
361     G4double XelPN  = csGGinstance->GetElasticIsotopeCrossSection(Proton, KineticEnergy, 0, 1);
362 
363     Xtotal   = ( NumberOfTargetProtons  * XtotPP + NumberOfTargetNeutrons * XtotPN )
364                / NumberOfTargetNucleons;
365     Xelastic = ( NumberOfTargetProtons  * XelPP  + NumberOfTargetNeutrons * XelPN )
366                / NumberOfTargetNucleons;
367     Xannihilation = 0.0;
368     Xtotal /= millibarn;
369     Xelastic /= millibarn;
370   };
371 
372   // Geometrical parameters
373   SetTotalCrossSection( Xtotal );
374   SetElasticCrossSection( Xelastic );
375   SetInelasticCrossSection( Xtotal - Xelastic );
376 
377   // Interactions with elastic and inelastic collisions
378   SetProbabilityOfElasticScatt( Xtotal, Xelastic );
379 
380   SetRadiusOfHNinteractions2( Xtotal/pi/10.0 );
381 
382   if ( ( Xtotal - Xelastic ) == 0.0 ) {
383     SetProbabilityOfAnnihilation( 0.0 );
384   } else {
385     SetProbabilityOfAnnihilation( Xannihilation / (Xtotal - Xelastic) );
386   }
387 
388   if(Xelastic > 0.0) {
389     SetSlope( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 );// Slope parameter of elastic scattering
390                                                       // (GeV/c)^(-2))
391     // Parameters of elastic scattering
392     // Gaussian parametrization of elastic scattering amplitude assumed
393     SetAvaragePt2ofElasticScattering( 1.0/( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 )*GeV*GeV );
394   } else {
395     SetSlope(1.0);
396     SetAvaragePt2ofElasticScattering( 0.0);
397   }
398   SetGamma0( GetSlope()*Xtotal/10.0/2.0/pi );
399 
400   G4double Xinel = Xtotal - Xelastic;
401 
402   #ifdef debugFTFparams
403   G4cout<< "Slope of hN elastic scattering" << GetSlope() << G4endl;
404   G4cout << "AvaragePt2ofElasticScattering " << GetAvaragePt2ofElasticScattering() << G4endl;
405   G4cout<<"Parameters of excitation for projectile "<<ProjectilePDGcode<< G4endl;
406   #endif
407 
408   if ( (ProjectilePDGcode == 2212) || (ProjectilePDGcode == 2112) ) {  // Projectile is proton or neutron
409     
410     const G4int indexTune = G4FTFTunings::Instance()->GetIndexTune( particle, KineticEnergy );
411     
412     // A process probability is parameterized as Prob = A_1*exp(-A_2*y) + A_3*exp(-A_4*y) + A_top
413     // y is a rapidity of a partcle in the target nucleus. Ymin is a minimal rapidity below it X=0
414 
415     //        Proc#   A1      B1            A2       B2   A3   Atop       Ymin
416     /* original hadr-string-diff-V10-03-07 (similar to 10.3.x) 
417     SetParams( 0,     13.71, 1.75,          -214.5, 4.25, 0.0, 0.5  ,     1.1 );  // Qexchange without Exc.
418     SetParams( 1,      25.0, 1.0,           -50.34, 1.5 , 0.0, 0.0  ,     1.4 );  // Qexchange with    Exc.
419     SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0  ,     0.93);  // Projectile diffraction
420     SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0  ,     0.93);  // Target diffraction
421     SetParams( 4,       1.0, 0.0 ,          -2.01 , 0.5 , 0.0, 0.0  ,     1.4 );  // Qexchange with Exc. Additional multiplier
422     */
423 
424     //         Proc#
425     SetParams( 0, fArrayParCollBaryonProj[indexTune].GetProc0A1(),
426             fArrayParCollBaryonProj[indexTune].GetProc0B1(),          
427       fArrayParCollBaryonProj[indexTune].GetProc0A2(),
428             fArrayParCollBaryonProj[indexTune].GetProc0B2(), 
429       fArrayParCollBaryonProj[indexTune].GetProc0A3(),
430             fArrayParCollBaryonProj[indexTune].GetProc0Atop(),     
431       fArrayParCollBaryonProj[indexTune].GetProc0Ymin() );   // Qexchange without Exc.
432     SetParams( 1, fArrayParCollBaryonProj[indexTune].GetProc1A1(),
433             fArrayParCollBaryonProj[indexTune].GetProc1B1(),          
434       fArrayParCollBaryonProj[indexTune].GetProc1A2(),
435             fArrayParCollBaryonProj[indexTune].GetProc1B2(), 
436       fArrayParCollBaryonProj[indexTune].GetProc1A3(),
437             fArrayParCollBaryonProj[indexTune].GetProc1Atop(),     
438       fArrayParCollBaryonProj[indexTune].GetProc1Ymin() );   // Qexchange with Exc.
439     if ( Xinel > 0.0 ) {
440       SetParams( 2, 6.0/Xinel, 0.0, -6.0/Xinel*16.28, 3.0, 0.0, 0.0, 0.93 );  // Projectile diffraction
441       SetParams( 3, 6.0/Xinel, 0.0, -6.0/Xinel*16.28, 3.0, 0.0, 0.0, 0.93 );  // Target diffraction
442 
443       SetParams( 4, fArrayParCollBaryonProj[indexTune].GetProc4A1(),
444         fArrayParCollBaryonProj[indexTune].GetProc4B1(),          
445         fArrayParCollBaryonProj[indexTune].GetProc4A2(),
446         fArrayParCollBaryonProj[indexTune].GetProc4B2(), 
447         fArrayParCollBaryonProj[indexTune].GetProc4A3(),
448         fArrayParCollBaryonProj[indexTune].GetProc4Atop(),     
449         fArrayParCollBaryonProj[indexTune].GetProc4Ymin() );  // Qexchange with Exc. Additional multiplier
450     } else {  // if Xinel=0., zero everything out (obviously)
451       SetParams( 2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
452       SetParams( 3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
453       SetParams( 4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
454     }
455 
456     if ( (AbsProjectileBaryonNumber > 10 || NumberOfTargetNucleons > 10) && !EnableDiffDissociationForBGreater10 ) {
457       // It is not decided what to do with diffraction dissociation in Had-Nucl and Nucl-Nucl interactions
458       // For the moment both ProjDiffDisso & TgtDiffDisso for A > 10 are set to false,
459       // so both projectile and target diffraction are turned OFF
460       if ( ! fArrayParCollBaryonProj[indexTune].IsProjDiffDissociation() )
461          SetParams( 2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -100.0 );  // Projectile diffraction
462       if ( ! fArrayParCollBaryonProj[indexTune].IsTgtDiffDissociation() )
463          SetParams( 3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -100.0 );  // Target diffraction
464     }
465 
466     SetDeltaProbAtQuarkExchange( fArrayParCollBaryonProj[indexTune].GetDeltaProbAtQuarkExchange() );
467 
468     if ( NumberOfTargetNucleons > 26 ) {
469       SetProbOfSameQuarkExchange( 1.0 );
470     } else {
471       SetProbOfSameQuarkExchange( fArrayParCollBaryonProj[indexTune].GetProbOfSameQuarkExchange() );
472     }
473 
474     SetProjMinDiffMass(    fArrayParCollBaryonProj[indexTune].GetProjMinDiffMass() );     // GeV 
475     SetProjMinNonDiffMass( fArrayParCollBaryonProj[indexTune].GetProjMinNonDiffMass() );  // GeV 
476 
477     SetTarMinDiffMass(     fArrayParCollBaryonProj[indexTune].GetTgtMinDiffMass() );      // GeV
478     SetTarMinNonDiffMass(  fArrayParCollBaryonProj[indexTune].GetTgtMinNonDiffMass() );   // GeV 
479 
480     SetAveragePt2(         fArrayParCollBaryonProj[indexTune].GetAveragePt2() );          // GeV^2       
481     SetProbLogDistrPrD(    fArrayParCollBaryonProj[indexTune].GetProbLogDistrPrD() ); 
482     SetProbLogDistr(       fArrayParCollBaryonProj[indexTune].GetProbLogDistr() ); 
483 
484   } else if ( ProjectilePDGcode == -2212  ||  ProjectilePDGcode == -2112 ) {  // Projectile is anti_proton or anti_neutron
485     
486     // Below, in the call to the G4FTFTunings::GetIndexTune method, we pass the proton
487     // as projectile, instead of the real one, because for switching on/off diffraction
488     // we assume the same treatment for anti_proton/anti_neutron as for proton/neutron,
489     // whereas all other parameters for anti_proton/anti_neutron are hardwired.
490     const G4int indexTune = G4FTFTunings::Instance()->GetIndexTune( G4Proton::Definition(), KineticEnergy );
491     
492     //        Proc#   A1      B1            A2       B2   A3   Atop       Ymin
493     SetParams( 0,      0.0 , 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,  1000.0  );  // Qexchange without Exc. 
494     SetParams( 1,      0.0 , 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,  1000.0  );  // Qexchange with    Exc.
495     if ( Xinel > 0.) {
496       SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 );  // Projectile diffraction
497       SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93 );  // Target diffraction
498       SetParams( 4,       1.0, 0.0 ,             0.0, 0.0 , 0.0, 0.0 , 0.93 );  // Qexchange with    Exc. Additional multiply
499     } else {
500       SetParams( 2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
501       SetParams( 3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
502       SetParams( 4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 );
503     }
504 
505     if ( AbsProjectileBaryonNumber > 10  ||  NumberOfTargetNucleons > 10 ) {
506       // It is not decided what to do with diffraction dissociation in Had-Nucl and Nucl-Nucl interactions
507       // For the moment both ProjDiffDisso & TgtDiffDisso are set to false,
508       // so both projectile and target diffraction are turned OFF
509       if ( ! fArrayParCollBaryonProj[indexTune].IsProjDiffDissociation() )
510          SetParams( 2,       0.0, 0.0 ,           0.0  , 0.0 , 0.0, 0.0   , -100.0  );  // Projectile diffraction
511       if ( ! fArrayParCollBaryonProj[indexTune].IsTgtDiffDissociation() )
512          SetParams( 3,       0.0, 0.0 ,           0.0  , 0.0 , 0.0, 0.0   , -100.0  );  // Target diffraction
513     }
514 
515     SetDeltaProbAtQuarkExchange( 0.0 );
516     SetProbOfSameQuarkExchange( 0.0 );
517     SetProjMinDiffMass( ProjectileMass + 0.22 );     // GeV 
518     SetProjMinNonDiffMass( ProjectileMass + 0.22 );  // GeV
519     SetTarMinDiffMass( TargetMass + 0.22 );          // GeV
520     SetTarMinNonDiffMass( TargetMass + 0.22 );       // GeV
521     SetAveragePt2( 0.3 );                            // GeV^2
522     SetProbLogDistrPrD( 0.55 );
523     SetProbLogDistr( 0.55 );
524             
525   } else if ( ProjectileabsPDGcode == 211  ||  ProjectilePDGcode ==  111 ) {  // Projectile is Pion
526     
527     const G4int indexTune = G4FTFTunings::Instance()->GetIndexTune( particle, KineticEnergy );
528  
529     //        Proc#   A1      B1            A2       B2      A3   Atop        Ymin
530     /* --> original code
531     SetParams( 0,  150.0,    1.8 ,       -247.3,     2.3,    0.,   1. ,       2.3 ); 
532     SetParams( 1,   5.77,    0.6 ,        -5.77,     0.8,    0.,   0. ,       0.0 );
533     SetParams( 2,   2.27,    0.5 ,     -98052.0,     4.0,    0.,   0. ,       3.0 ); 
534     SetParams( 3,    7.0,    0.9,        -85.28,     1.9,  0.08,   0. ,       2.2 );
535     SetParams( 4,    1.0,    0.0 ,       -11.02,     1.0,   0.0,   0. ,       2.4 );  // Qexchange with    Exc. Additional multiply
536     */
537     //         Proc#
538     SetParams( 0, fArrayParCollPionProj[indexTune].GetProc0A1(),
539             fArrayParCollPionProj[indexTune].GetProc0B1(),          
540       fArrayParCollPionProj[indexTune].GetProc0A2(),
541             fArrayParCollPionProj[indexTune].GetProc0B2(), 
542       fArrayParCollPionProj[indexTune].GetProc0A3(),
543             fArrayParCollPionProj[indexTune].GetProc0Atop(),     
544       fArrayParCollPionProj[indexTune].GetProc0Ymin() );   // Qexchange without Exc.
545     SetParams( 1, fArrayParCollPionProj[indexTune].GetProc1A1(),
546             fArrayParCollPionProj[indexTune].GetProc1B1(),          
547       fArrayParCollPionProj[indexTune].GetProc1A2(),
548             fArrayParCollPionProj[indexTune].GetProc1B2(), 
549       fArrayParCollPionProj[indexTune].GetProc1A3(),
550             fArrayParCollPionProj[indexTune].GetProc1Atop(),     
551       fArrayParCollPionProj[indexTune].GetProc1Ymin() );   // Qexchange with Exc.
552     SetParams( 2, fArrayParCollPionProj[indexTune].GetProc2A1(),
553             fArrayParCollPionProj[indexTune].GetProc2B1(),          
554       fArrayParCollPionProj[indexTune].GetProc2A2(),
555             fArrayParCollPionProj[indexTune].GetProc2B2(), 
556       fArrayParCollPionProj[indexTune].GetProc2A3(),
557             fArrayParCollPionProj[indexTune].GetProc2Atop(),     
558       fArrayParCollPionProj[indexTune].GetProc2Ymin() );   // Projectile diffraction
559     SetParams( 3, fArrayParCollPionProj[indexTune].GetProc3A1(),
560             fArrayParCollPionProj[indexTune].GetProc3B1(),          
561       fArrayParCollPionProj[indexTune].GetProc3A2(),
562             fArrayParCollPionProj[indexTune].GetProc3B2(), 
563       fArrayParCollPionProj[indexTune].GetProc3A3(),
564             fArrayParCollPionProj[indexTune].GetProc3Atop(),     
565       fArrayParCollPionProj[indexTune].GetProc3Ymin() );   // Target diffraction
566     SetParams( 4, fArrayParCollPionProj[indexTune].GetProc4A1(),
567             fArrayParCollPionProj[indexTune].GetProc4B1(),          
568       fArrayParCollPionProj[indexTune].GetProc4A2(),
569             fArrayParCollPionProj[indexTune].GetProc4B2(), 
570       fArrayParCollPionProj[indexTune].GetProc4A3(),
571             fArrayParCollPionProj[indexTune].GetProc4Atop(),     
572       fArrayParCollPionProj[indexTune].GetProc4Ymin() );   // Qexchange with Exc. Additional multiply
573 
574     // NOTE: how can it be |ProjectileBaryonNumber| > 10 if projectile is a pion ??? 
575     //
576     if ( AbsProjectileBaryonNumber > 10  ||  NumberOfTargetNucleons > 10 ) {
577        if ( ! fArrayParCollPionProj[indexTune].IsProjDiffDissociation() )
578           SetParams( 2,      0.0 , 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,  -100.0 );  // Projectile diffraction
579        if ( ! fArrayParCollPionProj[indexTune].IsTgtDiffDissociation() )
580           SetParams( 3,      0.0 , 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,  -100.0 );  // Target diffraction
581     }
582 
583     /* original code -->
584     SetDeltaProbAtQuarkExchange( 0.56 );
585     SetProjMinDiffMass( 1.0 );            // GeV
586     SetProjMinNonDiffMass( 1.0 );         // GeV 
587     SetTarMinDiffMass( 1.16 );            // GeV
588     SetTarMinNonDiffMass( 1.16 );         // GeV
589     SetAveragePt2( 0.3 );                 // GeV^2
590     SetProbLogDistrPrD( 0.55 );
591     SetProbLogDistr( 0.55 );
592     */
593 
594     // JVY update, Aug.8, 2018 --> Feb.14, 2019
595     //
596     SetDeltaProbAtQuarkExchange( fArrayParCollPionProj[indexTune].GetDeltaProbAtQuarkExchange() );
597     SetProjMinDiffMass( fArrayParCollPionProj[indexTune].GetProjMinDiffMass() );                  // GeV 
598     SetProjMinNonDiffMass( fArrayParCollPionProj[indexTune].GetProjMinNonDiffMass() );            // GeV 
599     SetTarMinDiffMass( fArrayParCollPionProj[indexTune].GetTgtMinDiffMass() );                    // GeV
600     SetTarMinNonDiffMass( fArrayParCollPionProj[indexTune].GetTgtMinNonDiffMass() );              // GeV 
601     SetAveragePt2( fArrayParCollPionProj[indexTune].GetAveragePt2() );                            // GeV^2       
602     SetProbLogDistrPrD( fArrayParCollPionProj[indexTune].GetProbLogDistrPrD() ); 
603     SetProbLogDistr( fArrayParCollPionProj[indexTune].GetProbLogDistr() ); 
604     
605     // ---> end update
606 
607   } else if ( ProjectileabsPDGcode == 321  ||  ProjectileabsPDGcode == 311  || 
608               ProjectilePDGcode == 130     ||  ProjectilePDGcode == 310 ) {  // Projectile is Kaon
609 
610     //        Proc#   A1      B1            A2       B2   A3   Atop       Ymin
611     SetParams( 0,     60.0 , 2.5 ,           0.0  , 0.0 , 0.0, 0.0  ,  -100.0  );  // Qexchange without Exc. 
612     SetParams( 1,      6.0 , 1.0 ,         -24.33 , 2.0 , 0.0, 0.0  ,     1.40 );  // Qexchange with    Exc.
613     SetParams( 2,      2.76, 1.2 ,         -22.5  , 2.7 ,0.04, 0.0  ,     1.40 );  // Projectile diffraction
614     SetParams( 3,      1.09, 0.5 ,          -8.88 , 2.  ,0.05, 0.0  ,     1.40 );  // Target diffraction
615     SetParams( 4,       1.0, 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,     0.93 );  // Qexchange with    Exc. Additional multiply
616     if ( AbsProjectileBaryonNumber > 10  ||  NumberOfTargetNucleons > 10 ) {
617       SetParams( 2,      0.0 , 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,  -100.0  );  // Projectile diffraction
618       SetParams( 3,      0.0 , 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,  -100.0  );  // Target diffraction
619     }
620 
621     SetDeltaProbAtQuarkExchange( 0.6 );
622     SetProjMinDiffMass( 0.7 );     // GeV 
623     SetProjMinNonDiffMass( 0.7 );  // GeV 
624     SetTarMinDiffMass( 1.16 );     // GeV
625     SetTarMinNonDiffMass( 1.16 );  // GeV
626     SetAveragePt2( 0.3 );          // GeV^2
627     SetProbLogDistrPrD( 0.55 );
628     SetProbLogDistr( 0.55 );
629 
630   } else {  // Projectile is not p, n, Pi0, Pi+, Pi-, K+, K-, K0 or their anti-particles
631 
632     if ( ProjectileabsPDGcode > 1000 ) {  // The projectile is a baryon as P or N
633       //        Proc#   A1      B1            A2       B2   A3   Atop       Ymin
634       SetParams( 0,     13.71, 1.75,          -30.69, 3.0 , 0.0, 1.0  ,     0.93 ); // Qexchange without Exc.
635       SetParams( 1,      25.0, 1.0,          -50.34,  1.5 , 0.0, 0.0  ,     1.4 );  // Qexchange with    Exc.
636       if ( Xinel > 0.) {
637         SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0  ,   0.93);  // Projectile diffraction
638         SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0  ,   0.93);  // Target diffraction
639         SetParams( 4,       1.0, 0.0 ,          -2.01 , 0.5 , 0.0, 0.0  ,   1.4 );  // Qexchange with    Exc. Additional multiply
640       } else {
641         SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0  ,     0.0);
642         SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0  ,     0.0);
643         SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0  ,     0.0);
644       }
645 
646     } else {  // The projectile is a meson as K+-0
647       //        Proc#   A1      B1            A2       B2   A3   Atop       Ymin
648       SetParams( 0,     60.0 , 2.5 ,           0.0  , 0.0 , 0.0, 0.0  ,  -100.0  );  // Qexchange without Exc. 
649       SetParams( 1,      6.0 , 1.0 ,         -24.33 , 2.0 , 0.0, 0.0  ,     1.40 );  // Qexchange with    Exc.
650       SetParams( 2,      2.76, 1.2 ,         -22.5  , 2.7 ,0.04, 0.0  ,     1.40 );  // Projectile diffraction
651       SetParams( 3,      1.09, 0.5 ,          -8.88 , 2.  ,0.05, 0.0  ,     1.40 );  // Target diffraction
652       SetParams( 4,       1.0, 0.0 ,           0.0  , 0.0 , 0.0, 0.0  ,     0.93 );  // Qexchange with    Exc. Additional multiply
653     }
654 
655     if ( AbsProjectileBaryonNumber > 10  ||  NumberOfTargetNucleons > 10 ) {
656       SetParams( 2,      0.0 , 0.0 ,            0.0 , 0.0 , 0.0, 0.0  ,  -100.0  );  // Projectile diffraction
657       SetParams( 3,      0.0 , 0.0 ,            0.0 , 0.0 , 0.0, 0.0  ,  -100.0  );  // Target diffraction
658     }
659 
660     SetDeltaProbAtQuarkExchange( 0.0 );
661     SetProbOfSameQuarkExchange( 0.0 );
662 
663     SetProjMinDiffMass(    GetMinMass(particle)/GeV );
664     SetProjMinNonDiffMass( GetMinMass(particle)/GeV );
665 
666     const G4ParticleDefinition* Neutron = G4Neutron::Neutron();
667     SetTarMinDiffMass(    GetMinMass(Neutron)/GeV );
668     SetTarMinNonDiffMass( GetMinMass(Neutron)/GeV );
669 
670     SetAveragePt2( 0.3 );  // GeV^2
671     SetProbLogDistrPrD( 0.55 );
672     SetProbLogDistr( 0.55 );
673 
674   }
675 
676   #ifdef debugFTFparams
677   G4cout<<"DeltaProbAtQuarkExchange "<< GetDeltaProbAtQuarkExchange() << G4endl;
678   G4cout<<"ProbOfSameQuarkExchange  "<< GetProbOfSameQuarkExchange()  << G4endl;
679   G4cout<<"ProjMinDiffMass          "<< GetProjMinDiffMass()/GeV <<" GeV"<< G4endl;
680   G4cout<<"ProjMinNonDiffMass       "<< GetProjMinNonDiffMass()  <<" GeV"<< G4endl;
681   G4cout<<"TarMinDiffMass           "<< GetTarMinDiffMass()      <<" GeV"<< G4endl;
682   G4cout<<"TarMinNonDiffMass        "<< GetTarMinNonDiffMass()   <<" GeV"<< G4endl;
683   G4cout<<"AveragePt2               "<< GetAveragePt2()          <<" GeV^2"<< G4endl;
684   G4cout<<"ProbLogDistrPrD          "<< GetProbLogDistrPrD() << G4endl;
685   G4cout<<"ProbLogDistrTrD          "<< GetProbLogDistr()    << G4endl;
686   #endif
687 
688   // Set parameters of nuclear destruction
689 
690   if ( ProjectileabsPDGcode < 1000 ) {  // Meson projectile
691 
692     const G4int indexTune = G4FTFTunings::Instance()->GetIndexTune( particle, KineticEnergy );
693 
694     SetMaxNumberOfCollisions( Plab, 2.0 );  //  3.0 )
695     //
696     // target destruction
697     //
698     /* original code --->
699     SetCofNuclearDestruction( 0.00481*G4double(NumberOfTargetNucleons)*
700             G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
701     
702     SetR2ofNuclearDestruction( 1.5*fermi*fermi );
703     SetDofNuclearDestruction( 0.3 );
704     SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
705                                          ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
706     SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
707     SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
708     */
709     double coeff = fArrayParCollMesonProj[indexTune].GetNuclearTgtDestructP1();
710     // 
711     // NOTE (JVY): Set this switch to false/true on line 138
712     //
713     if ( fArrayParCollMesonProj[indexTune].IsNuclearTgtDestructP1_ADEP() )       
714     {                                                             
715       coeff *= G4double(NumberOfTargetNucleons);                  
716     }                                                             
717     double exfactor = G4Exp( fArrayParCollMesonProj[indexTune].GetNuclearTgtDestructP2()
718                            * (Ylab-fArrayParCollMesonProj[indexTune].GetNuclearTgtDestructP3()) );
719     coeff *= exfactor;
720     coeff /= ( 1.+ exfactor );
721 
722     SetCofNuclearDestruction( coeff );
723 
724     SetR2ofNuclearDestruction( fArrayParCollMesonProj[indexTune].GetR2ofNuclearDestruct() );
725     SetDofNuclearDestruction( fArrayParCollMesonProj[indexTune].GetDofNuclearDestruct() );
726     coeff = fArrayParCollMesonProj[indexTune].GetPt2NuclearDestructP2();
727     exfactor = G4Exp( fArrayParCollMesonProj[indexTune].GetPt2NuclearDestructP3()
728                     * (Ylab-fArrayParCollMesonProj[indexTune].GetPt2NuclearDestructP4()) );
729     coeff *= exfactor;
730     coeff /= ( 1. + exfactor );
731     SetPt2ofNuclearDestruction( (fArrayParCollMesonProj[indexTune].GetPt2NuclearDestructP1()+coeff)*CLHEP::GeV*CLHEP::GeV ); 
732 
733     SetMaxPt2ofNuclearDestruction( fArrayParCollMesonProj[indexTune].GetMaxPt2ofNuclearDestruct() );
734     SetExcitationEnergyPerWoundedNucleon( fArrayParCollMesonProj[indexTune].GetExciEnergyPerWoundedNucleon() ); 
735 
736   } else if ( ProjectilePDGcode == -2212  ||  ProjectilePDGcode == -2112 ) {  // for anti-baryon projectile
737 
738     SetMaxNumberOfCollisions( Plab, 2.0 );
739 
740     SetCofNuclearDestruction( 0.00481*G4double(NumberOfTargetNucleons)*
741            G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
742     SetR2ofNuclearDestruction( 1.5*fermi*fermi );
743     SetDofNuclearDestruction( 0.3 );
744     SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
745                                          ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
746     SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
747     SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
748     if ( Plab < 2.0 ) {  // 2 GeV/c
749       // For slow anti-baryon we have to garanty putting on mass-shell
750       SetCofNuclearDestruction( 0.0 );
751       SetR2ofNuclearDestruction( 1.5*fermi*fermi ); // this is equivalent to setting a few line above
752                                                     // is it even necessary ?
753       SetDofNuclearDestruction( 0.01 );
754       SetPt2ofNuclearDestruction( 0.035*GeV*GeV );
755       SetMaxPt2ofNuclearDestruction( 0.04*GeV*GeV );
756     }
757 
758   } else {  // Projectile baryon assumed
759 
760     // Below, in the call to the G4FTFTunings::GetIndexTune method, we pass the proton
761     // as projectile, instead of the real one, because for the treatment of nuclear
762     // destruction, we assume for this category of hadron projectiles the same treatment
763     // as for "baryon".
764     const G4int indexTune = G4FTFTunings::Instance()->GetIndexTune( G4Proton::Definition(), KineticEnergy );
765     
766     // NOTE (JVY) FIXME !!! Will decide later how/if to make this one configurable...
767     //
768     SetMaxNumberOfCollisions( Plab, 2.0 );
769 
770     // projectile destruction - does NOT really matter for particle projectile, only for a nucleus projectile
771     //
772     double coeff = 0.;
773     coeff = fArrayParCollBaryonProj[indexTune].GetNuclearProjDestructP1();
774     // 
775     // NOTE (JVY): Set this switch to false/true on line 136
776     //
777     if ( fArrayParCollBaryonProj[indexTune].IsNuclearProjDestructP1_NBRNDEP() )   
778     {                                                             
779       coeff *= G4double(AbsProjectileBaryonNumber);               
780     }                                                             
781     double exfactor = G4Exp( fArrayParCollBaryonProj[indexTune].GetNuclearProjDestructP2()*
782            (Ylab-fArrayParCollBaryonProj[indexTune].GetNuclearProjDestructP3()) ); 
783     coeff *= exfactor;
784     coeff /= ( 1.+ exfactor );
785     SetCofNuclearDestructionPr( coeff );
786 
787     // target desctruction
788     //
789     coeff = fArrayParCollBaryonProj[indexTune].GetNuclearTgtDestructP1();
790     // 
791     // NOTE (JVY): Set this switch to false/true on line 138
792     //
793     if ( fArrayParCollBaryonProj[indexTune].IsNuclearTgtDestructP1_ADEP() )       
794     {                                                             
795       coeff *= G4double(NumberOfTargetNucleons);                  
796     }                                                             
797     exfactor = G4Exp( fArrayParCollBaryonProj[indexTune].GetNuclearTgtDestructP2()*
798           (Ylab-fArrayParCollBaryonProj[indexTune].GetNuclearTgtDestructP3()) );
799     coeff *= exfactor;
800     coeff /= ( 1.+ exfactor );
801     SetCofNuclearDestruction(   coeff );
802 
803     SetR2ofNuclearDestruction( fArrayParCollBaryonProj[indexTune].GetR2ofNuclearDestruct() );
804     SetDofNuclearDestruction( fArrayParCollBaryonProj[indexTune].GetDofNuclearDestruct() );
805 
806     coeff = fArrayParCollBaryonProj[indexTune].GetPt2NuclearDestructP2();
807     exfactor = G4Exp( fArrayParCollBaryonProj[indexTune].GetPt2NuclearDestructP3()*
808           (Ylab-fArrayParCollBaryonProj[indexTune].GetPt2NuclearDestructP4())  );
809     coeff *= exfactor;
810     coeff /= ( 1. + exfactor );
811     SetPt2ofNuclearDestruction( (fArrayParCollBaryonProj[indexTune].GetPt2NuclearDestructP1()+coeff)*CLHEP::GeV*CLHEP::GeV ); 
812 
813     SetMaxPt2ofNuclearDestruction( fArrayParCollBaryonProj[indexTune].GetMaxPt2ofNuclearDestruct() );
814     SetExcitationEnergyPerWoundedNucleon( fArrayParCollBaryonProj[indexTune].GetExciEnergyPerWoundedNucleon() ); 
815 
816   }
817 
818   #ifdef debugFTFparams
819   G4cout<<"CofNuclearDestructionPr           "<< GetCofNuclearDestructionPr() << G4endl;
820   G4cout<<"CofNuclearDestructionTr           "<< GetCofNuclearDestruction()   << G4endl;
821   G4cout<<"R2ofNuclearDestruction            "<< GetR2ofNuclearDestruction()/fermi/fermi  <<" fermi^2"<< G4endl;
822   G4cout<<"DofNuclearDestruction             "<< GetDofNuclearDestruction()  << G4endl;
823   G4cout<<"Pt2ofNuclearDestruction           "<< GetPt2ofNuclearDestruction()/GeV/GeV <<" GeV^2"<< G4endl;
824   G4cout<<"ExcitationEnergyPerWoundedNucleon "<< GetExcitationEnergyPerWoundedNucleon()   <<" MeV"<< G4endl;
825   #endif
826 
827   //SetCofNuclearDestruction( 0.47*G4Exp( 2.0*(Ylab - 2.5) )/( 1.0 + G4Exp( 2.0*(Ylab - 2.5) ) ) ); 
828   //SetPt2ofNuclearDestruction( ( 0.035 + 0.1*G4Exp( 4.0*(Ylab - 3.0) )/( 1.0 + G4Exp( 4.0*(Ylab - 3.0) ) ) )*GeV*GeV );
829 
830   //SetMagQuarkExchange( 120.0 );       // 210.0 PipP
831   //SetSlopeQuarkExchange( 2.0 );
832   //SetDeltaProbAtQuarkExchange( 0.6 );
833   //SetProjMinDiffMass( 0.7 );          // GeV 1.1
834   //SetProjMinNonDiffMass( 0.7 );       // GeV
835   //SetProbabilityOfProjDiff( 0.0);     // 0.85*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
836   //SetTarMinDiffMass( 1.1 );           // GeV
837   //SetTarMinNonDiffMass( 1.1 );        // GeV
838   //SetProbabilityOfTarDiff( 0.0 );     // 0.85*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
839 
840   //SetAveragePt2( 0.0 );               // GeV^2   0.3
841   //------------------------------------
842   //SetProbabilityOfElasticScatt( 1.0, 1.0);                            //(Xtotal, Xelastic);
843   //SetProbabilityOfProjDiff( 1.0*0.62*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.51 ) );  // 0->1
844   //SetProbabilityOfTarDiff( 4.0*0.62*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.51 ) );   // 2->4
845   //SetAveragePt2( 0.3 );                                               // (0.15)
846   //SetAvaragePt2ofElasticScattering( 0.0 );
847 
848   //SetMaxNumberOfCollisions( Plab, 6.0 ); //(4.0*(Plab + 0.01), Plab); // 6.0 );
849   //SetAveragePt2( 0.15 );
850   //SetCofNuclearDestruction(-1.);//( 0.75 );                           // (0.25)             
851   //SetExcitationEnergyPerWoundedNucleon(0.);//( 30.0*MeV );            // (75.0*MeV) 
852   //SetDofNuclearDestruction(0.);//( 0.2 ); //0.4                       // 0.3 0.5
853 
854   /*
855   SetAveragePt2(0.3);
856   SetCofNuclearDestructionPr(0.);
857   SetCofNuclearDestruction(0.);             //( 0.5 ); (0.25)             
858   SetExcitationEnergyPerWoundedNucleon(0.); // 30.0*MeV; (75.0*MeV) 
859   SetDofNuclearDestruction(0.);             // 0.2; 0.4; 0.3; 0.5
860   SetPt2ofNuclearDestruction(0.);           //(2.*0.075*GeV*GeV); ( 0.3*GeV*GeV ); (0.168*GeV*GeV) 
861   */
862 
863   //SetExcitationEnergyPerWoundedNucleon(0.001);
864   //SetPt2Kink( 0.0*GeV*GeV );
865 
866   //SetRadiusOfHNinteractions2( Xtotal/pi/10.0 /2.);
867   //SetRadiusOfHNinteractions2( (Xtotal - Xelastic)/pi/10.0 );
868   //SetProbabilityOfElasticScatt( 1.0, 0.0);
869   /*
870   G4cout << "Pt2 " << GetAveragePt2()<<" "<<GetAveragePt2()/GeV/GeV<<G4endl;
871   G4cout << "Cnd " << GetCofNuclearDestruction() << G4endl;
872   G4cout << "Dnd " << GetDofNuclearDestruction() << G4endl;
873   G4cout << "Pt2 " << GetPt2ofNuclearDestruction()/GeV/GeV << G4endl;
874   */
875 
876 }
877 
878 //============================================================================
879 
880 G4double G4FTFParameters::GetMinMass( const G4ParticleDefinition* aParticle ) {
881    // The code is used for estimating the minimal string mass produced in diffraction dissociation.
882    // The indices used for minMassQDiQStr must be between 1 and 5, corresponding to the 5 considered
883    // quarks: d, u, s, c and b; enforcing this explicitly avoids compilation errors.
884    G4double EstimatedMass = 0.0;
885    G4int partID = std::abs(aParticle->GetPDGEncoding());
886    G4int Qleft  = std::max(  partID/100,     1 );
887    G4int Qright = std::max( (partID/ 10)%10, 1 );
888    if        ( Qleft < 6  &&  Qright < 6 ) {  // Q-Qbar string
889      EstimatedMass = StringMass->minMassQQbarStr[Qleft-1][Qright-1];
890    } else if ( Qleft < 6  &&  Qright > 6 ) {  // Q - DiQ string
891      G4int q1 = std::max( std::min( Qright/10, 5 ), 1 );
892      G4int q2 = std::max( std::min( Qright%10, 5 ), 1 );
893      EstimatedMass = StringMass->minMassQDiQStr[Qleft-1][q1-1][q2-1];
894    } else if ( Qleft > 6  &&  Qright < 6 ) {  // DiQ - Q string
895      G4int q1 = std::max( std::min( Qleft/10, 5 ), 1 );
896      G4int q2 = std::max( std::min( Qleft%10, 5 ), 1 );
897      EstimatedMass = StringMass->minMassQDiQStr[Qright-1][q1-1][q2-1];
898    }
899    return EstimatedMass;
900 }
901 
902 //============================================================================
903 
904 G4double G4FTFParameters::GetProcProb( const G4int ProcN, const G4double y ) {
905   G4double Prob( 0.0 );
906   if ( y < ProcParams[ProcN][6] ) {
907     Prob = ProcParams[ProcN][5]; 
908     if (Prob < 0.) Prob=0.;
909     return Prob;
910   }
911   Prob = ProcParams[ProcN][0] * G4Exp( -ProcParams[ProcN][1]*y ) +
912          ProcParams[ProcN][2] * G4Exp( -ProcParams[ProcN][3]*y ) +
913          ProcParams[ProcN][4];
914   if (Prob < 0.) Prob=0.;
915   return Prob;
916 }
917 
918 //============================================================================
919 
920 G4FTFParameters::~G4FTFParameters() {
921   if ( StringMass ) delete StringMass;
922 }
923 
924 //============================================================================
925 
926 void G4FTFParameters::Reset()
927 {
928   FTFhNcmsEnergy  = 0.0; 
929   FTFXtotal = 0.0;
930   FTFXelastic = 0.0;
931   FTFXinelastic = 0.0; 
932   FTFXannihilation = 0.0;
933   ProbabilityOfAnnihilation = 0.0; 
934   ProbabilityOfElasticScatt  = 0.0;
935   RadiusOfHNinteractions2 = 0.0; 
936   FTFSlope = 0.0;  
937   AvaragePt2ofElasticScattering = 0.0; 
938   FTFGamma0 = 0.0;
939   DeltaProbAtQuarkExchange = 0.0; 
940   ProbOfSameQuarkExchange = 0.0; 
941   ProjMinDiffMass = 0.0; 
942   ProjMinNonDiffMass = 0.0; 
943   ProbLogDistrPrD = 0.0;
944   TarMinDiffMass = 0.0; 
945   TarMinNonDiffMass = 0.0;
946   AveragePt2 = 0.0; 
947   ProbLogDistr  = 0.0;
948   Pt2kink = 0.0;
949   MaxNumberOfCollisions = 0.0;  
950   ProbOfInelInteraction = 0.0; 
951   CofNuclearDestructionPr = 0.0; 
952   CofNuclearDestruction = 0.0;
953   R2ofNuclearDestruction  = 0.0; 
954   ExcitationEnergyPerWoundedNucleon  = 0.0;
955   DofNuclearDestruction  = 0.0; 
956   Pt2ofNuclearDestruction = 0.0; 
957   MaxPt2ofNuclearDestruction = 0.0; 
958 
959   for ( G4int i = 0; i < 4; i++ ) {
960     for ( G4int j = 0; j < 7; j++ ) {
961       ProcParams[i][j] = 0.0;
962     }
963   }
964 
965   return;
966 }
967 
968