Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FTFParameters.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id$
28 // GEANT4 tag $Name: $
29 //
30 
31 #include <utility>
32 
33 #include "G4FTFParameters.hh"
34 
35 #include "G4ios.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 
39 #include "G4ParticleDefinition.hh" // 31 May 2011
40 
41 #include "G4Proton.hh" // 31 May 2011
42 #include "G4Neutron.hh" // 31 May 2011
43 
44 #include "G4PionPlus.hh" // 31 May 2011
45 #include "G4PionMinus.hh" // 31 May 2011
46 #include "G4KaonPlus.hh" // 31 May 2011
47 #include "G4KaonMinus.hh" // 31 May 2011
48 
49 G4FTFParameters::G4FTFParameters() :
50  //A.R. 14-Aug-2012 : Coverity fix.
51  FTFhNcmsEnergy(0.0),
52  FTFxsManager(0),
53  FTFXtotal(0.0), FTFXelastic(0.0), FTFXinelastic(0.0), FTFXannihilation(0.0),
54  ProbabilityOfAnnihilation(0.0), ProbabilityOfElasticScatt(0.0),
55  RadiusOfHNinteractions2(0.0), FTFSlope(0.0),
56  AvaragePt2ofElasticScattering(0.0), FTFGamma0(0.0),
57  MagQuarkExchange(0.0), SlopeQuarkExchange(0.0), DeltaProbAtQuarkExchange(0.0),
58  ProbOfSameQuarkExchange(0.0), ProjMinDiffMass(0.0), ProjMinNonDiffMass(0.0),
59  ProbabilityOfProjDiff(0.0), TarMinDiffMass(0.0), TarMinNonDiffMass(0.0),
60  ProbabilityOfTarDiff(0.0), AveragePt2(0.0), ProbLogDistr(0.0),
61  Pt2kink(0.0),
62  MaxNumberOfCollisions(0.0), ProbOfInelInteraction(0.0), CofNuclearDestruction(0.0),
63  R2ofNuclearDestruction(0.0), ExcitationEnergyPerWoundedNucleon(0.0),
64  DofNuclearDestruction(0.0), Pt2ofNuclearDestruction(0.0), MaxPt2ofNuclearDestruction(0.0)
65 {}
66 
67 
68 G4FTFParameters::~G4FTFParameters()
69 {}
70 //**********************************************************************************************
71 G4FTFParameters::G4FTFParameters(const G4ParticleDefinition * particle,
72  G4int theA,
73  G4int theZ,
74  G4double PlabPerParticle)
75 {
76 
77  //A.R. 25-Jul-2012 Coverity fix.
78  FTFXannihilation = 0.0;
79  FTFhNcmsEnergy = 0.0;
80  ProbOfSameQuarkExchange = 0.0;
81 
82  G4int ProjectilePDGcode = particle->GetPDGEncoding();
83  G4int ProjectileabsPDGcode = std::abs(ProjectilePDGcode);
84  G4double ProjectileMass = particle->GetPDGMass();
85  G4double ProjectileMass2 =ProjectileMass*ProjectileMass;
86 
87  G4int ProjectileBaryonNumber(0), AbsProjectileBaryonNumber(0);
88  G4int AbsProjectileCharge(0);
89  G4bool ProjectileIsNucleus=false;
90 
91  if(std::abs(particle->GetBaryonNumber()) > 1)
92  { // The projectile is a nucleus
93  ProjectileIsNucleus =true;
94  ProjectileBaryonNumber =particle->GetBaryonNumber();
95  AbsProjectileBaryonNumber=std::abs(ProjectileBaryonNumber);
96  AbsProjectileCharge =(G4int) particle->GetPDGCharge();
97 
98  if(ProjectileBaryonNumber > 1)
99  { ProjectilePDGcode= 2212; ProjectileabsPDGcode=2212;} // Proton
100  else { ProjectilePDGcode=-2212; ProjectileabsPDGcode=2212;} // Anti-Proton
101 
102  ProjectileMass =G4Proton::Proton()->GetPDGMass();
103  ProjectileMass2=sqr(ProjectileMass);
104  }
105 
106  G4double TargetMass = G4Proton::Proton()->GetPDGMass();
107  G4double TargetMass2 = TargetMass*TargetMass;
108 
109  G4double Plab = PlabPerParticle;
110  G4double Elab = std::sqrt(Plab*Plab+ProjectileMass2);
111  G4double KineticEnergy = Elab-ProjectileMass; // 31 May 2011
112 
113  G4double S=ProjectileMass2 + TargetMass2 + 2.*TargetMass*Elab;
114 
115 //G4cout<<"Proj Plab "<<ProjectilePDGcode<<" "<<Plab<<G4endl;
116 //G4cout<<"Mass KinE "<<ProjectileMass<<" "<<KineticEnergy<<G4endl;
117 //G4cout<<" A Z "<<theA<<" "<<theZ<<G4endl;
118 
119  G4double Ylab,Xtotal,Xelastic,Xannihilation;
120  G4int NumberOfTargetNucleons;
121 
122  Ylab=0.5*std::log((Elab+Plab)/(Elab-Plab));
123 
124  G4double ECMSsqr=S/GeV/GeV;
125  G4double SqrtS =std::sqrt(S)/GeV;
126 //G4cout<<"Sqrt(s) "<<SqrtS<<G4endl;
127 
128  TargetMass /=GeV; TargetMass2 /=(GeV*GeV);
129  ProjectileMass /=GeV; ProjectileMass2 /=(GeV*GeV);
130 
131  static G4ChipsComponentXS* _instance = new G4ChipsComponentXS(); // Witek Pokorski
132  FTFxsManager = _instance;
133 
134  Plab/=GeV;
135 // G4double LogPlab = std::log( Plab );
136 // G4double sqrLogPlab = LogPlab * LogPlab;
137 
138  G4int NumberOfTargetProtons = theZ;
139  G4int NumberOfTargetNeutrons = theA-theZ;
140 
141  NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
142 
143  if( (ProjectilePDGcode == 2212) ||
144  (ProjectilePDGcode == 2112) ) //------Projectile is nucleon --------
145  {
146  G4double XtotPP = FTFxsManager->
147  GetTotalElementCrossSection( particle,KineticEnergy,1,0);
148  G4ParticleDefinition* Neutron=G4Neutron::Neutron();
149  G4double XtotPN = FTFxsManager->
150  GetTotalElementCrossSection( Neutron,KineticEnergy,1,0);
151 
152 
153  G4double XelPP = FTFxsManager->
154  GetElasticElementCrossSection(particle,KineticEnergy,1,0);
155  G4double XelPN = FTFxsManager->
156  GetElasticElementCrossSection( Neutron,KineticEnergy,1,0);
157 //G4cout<<"Xs "<<XtotPP/millibarn<<" "<<XelPP/millibarn<<G4endl;
158 //G4cout<<"Xs "<<XtotPN/millibarn<<" "<<XelPN/millibarn<<G4endl;
159  if(!ProjectileIsNucleus)
160  { // Projectile is hadron
161  Xtotal = ( NumberOfTargetProtons * XtotPP +
162  NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
163  Xelastic = ( NumberOfTargetProtons * XelPP +
164  NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
165  } else
166  { // Projectile is a nucleus
167  Xtotal = (
168  AbsProjectileCharge *NumberOfTargetProtons *XtotPP +
169  (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetNeutrons*XtotPP +
170  ( AbsProjectileCharge *NumberOfTargetNeutrons +
171  (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetProtons)*XtotPN
172  )/(AbsProjectileBaryonNumber*NumberOfTargetNucleons);
173 
174  Xelastic= (
175  AbsProjectileCharge *NumberOfTargetProtons *XelPP +
176  (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetNeutrons*XelPP +
177  ( AbsProjectileCharge *NumberOfTargetNeutrons +
178  (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetProtons)*XelPN
179  )/(AbsProjectileBaryonNumber*NumberOfTargetNucleons);
180  }
181 
182  Xannihilation = 0.;
183  Xtotal/=millibarn;
184  Xelastic/=millibarn;
185  }
186  else if( ProjectilePDGcode < -1000 ) //------Projectile is anti_baryon --------
187  {
188 
189  G4double X_a(0.), X_b(0.), X_c(0.), X_d(0.);
190  G4double MesonProdThreshold=ProjectileMass+TargetMass+(2.*0.14+0.016); // 2 Mpi +DeltaE;
191 
192  if(PlabPerParticle < 40.*MeV)
193  { // Low energy limits. Projectile at rest.
194  Xtotal= 1512.9; // mb
195  Xelastic= 473.2; // mb
196  X_a= 625.1; // mb
197  X_b= 9.780; // mb
198  X_c= 49.989; // mb
199  X_d= 6.614; // mb
200  }
201  else
202  { // Total and elastic cross section of PbarP interactions a'la Arkhipov
203  G4double LogS=std::log(ECMSsqr/33.0625);
204  G4double Xasmpt=36.04+0.304*LogS*LogS; // mb
205 
206  LogS=std::log(SqrtS/20.74);
207  G4double Basmpt=11.92+0.3036*LogS*LogS; // GeV^(-2)
208  G4double R0=std::sqrt(0.40874044*Xasmpt-Basmpt); // GeV^(-1)
209 
210  G4double FlowF=SqrtS/
211  std::sqrt(ECMSsqr*ECMSsqr+ProjectileMass2*ProjectileMass2+TargetMass2*TargetMass2-
212  2.*ECMSsqr*ProjectileMass2 -2.*ECMSsqr*TargetMass2 -2.*ProjectileMass2*TargetMass2);
213 
214  Xtotal=Xasmpt*(1.+13.55*FlowF/R0/R0/R0*
215  (1.-4.47/SqrtS+12.38/ECMSsqr-12.43/SqrtS/ECMSsqr)); // mb
216 
217  Xasmpt=4.4+0.101*LogS*LogS; // mb
218  Xelastic=Xasmpt*(1.+59.27*FlowF/R0/R0/R0*
219  (1.-6.95/SqrtS+23.54/ECMSsqr-25.34/SqrtS/ECMSsqr)); // mb
220 //G4cout<<"Param Xtotal Xelastic "<<Xtotal<<" "<<Xelastic<<G4endl;
221 //G4cout<<"FlowF "<<FlowF<<" SqrtS "<<SqrtS<<G4endl;
222 //G4cout<<"Param Xelastic-NaN "<<Xelastic<<" "<<1.5*16.654/pow(ECMSsqr/2.176/2.176,2.2)<<" "<<ECMSsqr<<G4endl;
223  X_a=25.*FlowF; // mb, 3-shirts diagram
224 
225  if(SqrtS < MesonProdThreshold)
226  {
227  X_b=3.13+140.*std::pow(MesonProdThreshold-SqrtS,2.5);// mb anti-quark-quark annihilation
228  Xelastic-=3.*X_b; // Xel-X(PbarP->NNbar)
229  } else
230  {
231  X_b=6.8/SqrtS; // mb anti-quark-quark annihilation
232  Xelastic-=3.*X_b; // Xel-X(PbarP->NNbar)
233  }
234 
235  X_c=2.*FlowF*sqr(ProjectileMass+TargetMass)/ECMSsqr; // mb rearrangement
236 
237 //G4cout<<"Old new Xa "<<35.*FlowF<<" "<<25.*FlowF<<G4endl;
238 
239  X_d=23.3/ECMSsqr; // mb anti-quark-quark string creation
240  }
241 //---------------------------------------------------------------
242 //G4cout<<"Param Xtotal Xelastic "<<Xtotal<<" "<<Xelastic<<G4endl;
243 //G4cout<<"Para a b c d "<<X_a<<" "<<X_b<<" "<<X_c<<" "<<X_d<<G4endl;
244 //G4cout<<"Para a b c d "<<X_a<<" "<<5.*X_b<<" "<<5.*X_c<<" "<<6.*X_d<<G4endl;
245  G4double Xann_on_P(0.), Xann_on_N(0.);
246 
247  if(ProjectilePDGcode == -2212) // Pbar+P/N
248  {Xann_on_P=X_a + X_b*5. + X_c*5. + X_d*6.; Xann_on_N=X_a + X_b*4. + X_c*4. + X_d*4.;}
249  else if(ProjectilePDGcode == -2112) // NeutrBar+P/N
250  {Xann_on_P=X_a + X_b*4. + X_c*4. + X_d*4.; Xann_on_N=X_a + X_b*5. + X_c*5. + X_d*6.;}
251  else if(ProjectilePDGcode == -3122) // LambdaBar+P/N
252  {Xann_on_P=X_a + X_b*3. + X_c*3. + X_d*2.; Xann_on_N=X_a + X_b*3. + X_c*3. + X_d*2.;}
253  else if(ProjectilePDGcode == -3112) // Sigma-Bar+P/N
254  {Xann_on_P=X_a + X_b*2. + X_c*2. + X_d*0.; Xann_on_N=X_a + X_b*4. + X_c*4. + X_d*2.;}
255  else if(ProjectilePDGcode == -3212) // Sigma0Bar+P/N
256  {Xann_on_P=X_a + X_b*3. + X_c*3. + X_d*2.; Xann_on_N=X_a + X_b*3. + X_c*3. + X_d*2.;}
257  else if(ProjectilePDGcode == -3222) // Sigma+Bar+P/N
258  {Xann_on_P=X_a + X_b*4. + X_c*4. + X_d*2.; Xann_on_N=X_a + X_b*2. + X_c*2. + X_d*0.;}
259  else if(ProjectilePDGcode == -3312) // Xi-Bar+P/N
260  {Xann_on_P=X_a + X_b*1. + X_c*1. + X_d*0.; Xann_on_N=X_a + X_b*2. + X_c*2. + X_d*0.;}
261  else if(ProjectilePDGcode == -3322) // Xi0Bar+P/N
262  {Xann_on_P=X_a + X_b*2. + X_c*2. + X_d*0.; Xann_on_N=X_a + X_b*1. + X_c*1. + X_d*0.;}
263  else if(ProjectilePDGcode == -3334) // Omega-Bar+P/N
264  {Xann_on_P=X_a + X_b*0. + X_c*0. + X_d*0.; Xann_on_N=X_a + X_b*0. + X_c*0. + X_d*0.;}
265  else {G4cout<<"Unknown anti-baryon for FTF annihilation"<<G4endl;}
266 //---------------------------------------------------------------
267 
268 //G4cout<<"Sum "<<Xann_on_P<<G4endl;
269 
270  if(!ProjectileIsNucleus)
271  { // Projectile is anti-baryon
272  Xannihilation = ( NumberOfTargetProtons * Xann_on_P +
273  NumberOfTargetNeutrons * Xann_on_N ) / NumberOfTargetNucleons;
274  } else
275  { // Projectile is a nucleus
276  Xannihilation=(
277  ( AbsProjectileCharge *NumberOfTargetProtons+
278  (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetNeutrons )*Xann_on_P +
279  ( AbsProjectileCharge *NumberOfTargetNeutrons+
280  (AbsProjectileBaryonNumber-AbsProjectileCharge)*NumberOfTargetProtons )*Xann_on_N
281  )/(AbsProjectileBaryonNumber*NumberOfTargetNucleons);
282  }
283 
284  G4double Xftf=0.;
285  MesonProdThreshold=ProjectileMass+TargetMass+(0.14+0.08); // Mpi +DeltaE
286  if(SqrtS > MesonProdThreshold) {Xftf=36.*(1.-MesonProdThreshold/SqrtS);}
287 
288  Xtotal = Xelastic + Xannihilation + Xftf;
289 /*
290 G4cout<<"Plab Xtotal, Xelastic Xinel Xftf "<<Plab<<" "<<Xtotal<<" "<<Xelastic<<" "<<Xtotal-Xelastic<<" "<<Xtotal-Xelastic-Xannihilation<<G4endl;
291 G4cout<<"Plab Xelastic/Xtotal, Xann/Xin "<<Plab<<" "<<Xelastic/Xtotal<<" "<<Xannihilation/(Xtotal-Xelastic)<<G4endl;
292 //G4int Uzhi; G4cin>>Uzhi;
293 */
294 //---------------------------------------------------------------
295  }
296  else if( ProjectilePDGcode == 211 ) //------Projectile is PionPlus -------
297  {
298  G4double XtotPiP = FTFxsManager->
299  GetTotalElementCrossSection( particle,KineticEnergy,1,0);
300  G4ParticleDefinition* PionMinus=G4PionMinus::PionMinus();
301  G4double XtotPiN = FTFxsManager->
302  GetTotalElementCrossSection( PionMinus,KineticEnergy,1,0);
303 
304  G4double XelPiP = FTFxsManager->
305  GetElasticElementCrossSection(particle,KineticEnergy,1,0);
306  G4double XelPiN = FTFxsManager->
307  GetElasticElementCrossSection( PionMinus,KineticEnergy,1,0);
308 
309  Xtotal = ( NumberOfTargetProtons * XtotPiP +
310  NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
311  Xelastic = ( NumberOfTargetProtons * XelPiP +
312  NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
313  Xannihilation = 0.;
314  Xtotal/=millibarn;
315  Xelastic/=millibarn;
316  }
317  else if( ProjectilePDGcode == -211 ) //------Projectile is PionMinus -------
318  {
319  G4double XtotPiP = FTFxsManager->
320  GetTotalElementCrossSection( particle,KineticEnergy,1,0);
321  G4ParticleDefinition* PionPlus=G4PionPlus::PionPlus();
322  G4double XtotPiN = FTFxsManager->
323  GetTotalElementCrossSection( PionPlus,KineticEnergy,1,0);
324 
325  G4double XelPiP = FTFxsManager->
326  GetElasticElementCrossSection(particle,KineticEnergy,1,0);
327  G4double XelPiN = FTFxsManager->
328  GetElasticElementCrossSection( PionPlus,KineticEnergy,1,0);
329 
330  Xtotal = ( NumberOfTargetProtons * XtotPiP +
331  NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
332  Xelastic = ( NumberOfTargetProtons * XelPiP +
333  NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
334  Xannihilation = 0.;
335  Xtotal/=millibarn;
336  Xelastic/=millibarn;
337  }
338 
339  else if( ProjectilePDGcode == 111 ) //------Projectile is PionZero -------
340  {
341  G4ParticleDefinition* PionPlus=G4PionPlus::PionPlus();
342  G4double XtotPipP= FTFxsManager->
343  GetTotalElementCrossSection( PionPlus,KineticEnergy,1,0);
344 
345  G4ParticleDefinition* PionMinus=G4PionMinus::PionMinus();
346  G4double XtotPimP= FTFxsManager->
347  GetTotalElementCrossSection( PionMinus,KineticEnergy,1,0);
348 
349  G4double XelPipP = FTFxsManager->
350  GetElasticElementCrossSection( PionPlus,KineticEnergy,1,0);
351  G4double XelPimP = FTFxsManager->
352  GetElasticElementCrossSection( PionMinus,KineticEnergy,1,0);
353 
354  G4double XtotPiP= (XtotPipP + XtotPimP)/2.;
355  G4double XtotPiN=XtotPiP;
356  G4double XelPiP = (XelPipP + XelPimP )/2.;
357  G4double XelPiN = XelPiP;
358 
359  Xtotal = ( NumberOfTargetProtons * XtotPiP +
360  NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
361  Xelastic = ( NumberOfTargetProtons * XelPiP +
362  NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
363  Xannihilation = 0.;
364 
365  Xtotal/=millibarn;
366  Xelastic/=millibarn;
367  }
368  else if( ProjectilePDGcode == 321 ) //------Projectile is KaonPlus -------
369  {
370  G4double XtotKP = FTFxsManager->
371  GetTotalElementCrossSection( particle,KineticEnergy,1,0);
372 
373  G4ParticleDefinition* KaonMinus=G4KaonMinus::KaonMinus();
374  G4double XtotKN = FTFxsManager->
375  GetTotalElementCrossSection( KaonMinus,KineticEnergy,1,0);
376 
377  G4double XelKP = FTFxsManager->
378  GetElasticElementCrossSection(particle,KineticEnergy,1,0);
379  G4double XelKN = FTFxsManager->
380  GetElasticElementCrossSection( KaonMinus,KineticEnergy,1,0);
381 
382  Xtotal = ( NumberOfTargetProtons * XtotKP +
383  NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
384  Xelastic = ( NumberOfTargetProtons * XelKP +
385  NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
386  Xannihilation = 0.;
387 
388  Xtotal/=millibarn;
389  Xelastic/=millibarn;
390  }
391  else if( ProjectilePDGcode ==-321 ) //------Projectile is KaonMinus ------
392  {
393  G4double XtotKP = FTFxsManager->
394  GetTotalElementCrossSection( particle,KineticEnergy,1,0);
395 
396  G4ParticleDefinition* KaonPlus=G4KaonPlus::KaonPlus();
397  G4double XtotKN = FTFxsManager->
398  GetTotalElementCrossSection( KaonPlus,KineticEnergy,1,0);
399 
400  G4double XelKP = FTFxsManager->
401  GetElasticElementCrossSection(particle,KineticEnergy,1,0);
402 
403  G4double XelKN = FTFxsManager->
404  GetElasticElementCrossSection(KaonPlus,KineticEnergy,1,0);
405 
406  Xtotal = ( NumberOfTargetProtons * XtotKP +
407  NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
408  Xelastic = ( NumberOfTargetProtons * XelKP +
409  NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
410  Xannihilation = 0.;
411 
412  Xtotal/=millibarn;
413  Xelastic/=millibarn;
414  }
415  else if((ProjectilePDGcode == 311) ||
416  (ProjectilePDGcode == 130) ||
417  (ProjectilePDGcode == 310)) //Projectile is KaonZero
418  {
419  G4ParticleDefinition* KaonPlus=G4KaonPlus::KaonPlus();
420  G4double XtotKpP= FTFxsManager->
421  GetTotalElementCrossSection( KaonPlus,KineticEnergy,1,0);
422 
423  G4ParticleDefinition* KaonMinus=G4KaonMinus::KaonMinus();
424  G4double XtotKmP= FTFxsManager->
425  GetTotalElementCrossSection( KaonMinus,KineticEnergy,1,0);
426 
427  G4double XelKpP = FTFxsManager->
428  GetElasticElementCrossSection( KaonPlus,KineticEnergy,1,0);
429  G4double XelKmP = FTFxsManager->
430  GetElasticElementCrossSection( KaonMinus,KineticEnergy,1,0);
431 
432  G4double XtotKP=(XtotKpP+XtotKmP)/2.;
433  G4double XtotKN=XtotKP;
434  G4double XelKP =(XelKpP +XelKmP )/2.;
435  G4double XelKN =XelKP;
436 
437  Xtotal = ( NumberOfTargetProtons * XtotKP +
438  NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
439  Xelastic = ( NumberOfTargetProtons * XelKP +
440  NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
441  Xannihilation = 0.;
442 
443  Xtotal/=millibarn;
444  Xelastic/=millibarn;
445  }
446  else //------Projectile is undefined, Nucleon assumed
447  {
448  G4ParticleDefinition* Proton=G4Proton::Proton();
449  G4double XtotPP = FTFxsManager->
450  GetTotalElementCrossSection( Proton,KineticEnergy,1,0);
451 
452  G4ParticleDefinition* Neutron=G4Neutron::Neutron();
453  G4double XtotPN = FTFxsManager->
454  GetTotalElementCrossSection( Neutron,KineticEnergy,1,0);
455 
456  G4double XelPP = FTFxsManager->
457  GetElasticElementCrossSection(Proton,KineticEnergy,1,0);
458  G4double XelPN = FTFxsManager->
459  GetElasticElementCrossSection( Neutron,KineticEnergy,1,0);
460 
461  Xtotal = ( NumberOfTargetProtons * XtotPP +
462  NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
463  Xelastic = ( NumberOfTargetProtons * XelPP +
464  NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
465  Xannihilation = 0.;
466 
467  Xtotal/=millibarn;
468  Xelastic/=millibarn;
469  };
470 
471 //----------- Geometrical parameters ------------------------------------------------
472  SetTotalCrossSection(Xtotal);
473  SetElastisCrossSection(Xelastic);
474  SetInelasticCrossSection(Xtotal-Xelastic);
475 
476 /*
477 G4cout<<"Plab Xtotal, Xelastic Xinel Xftf "<<Plab<<" "<<Xtotal<<" "<<Xelastic<<" "<<Xtotal-Xelastic<<" "<<Xtotal-Xelastic-Xannihilation<<G4endl;
478 if(Xtotal-Xelastic != 0.)
479 {
480  G4cout<<"Plab Xelastic/Xtotal, Xann/Xin "<<Plab<<" "<<Xelastic/Xtotal<<" "<<Xannihilation/
481  (Xtotal-Xelastic)<<G4endl;
482 } else
483 {
484  G4cout<<"Plab Xelastic/Xtotal, Xann "<<Plab<<" "<<Xelastic/Xtotal<<" "<<
485  Xannihilation<<G4endl;
486 }
487 //G4int Uzhi; G4cin>>Uzhi;
488 */
489 // // Interactions with elastic and inelastic collisions
490  SetProbabilityOfElasticScatt(Xtotal, Xelastic);
491  SetRadiusOfHNinteractions2(Xtotal/pi/10.);
492 
493  if(Xtotal-Xelastic == 0.)
494  {
495  SetProbabilityOfAnnihilation(0.);
496  } else
497  {SetProbabilityOfAnnihilation(Xannihilation/(Xtotal-Xelastic));}
498 //
499 //SetProbabilityOfElasticScatt(Xtotal, 0.);
500 // //==== No elastic scattering ============================
501 // SetProbabilityOfElasticScatt(Xtotal, 0.);
502 // SetRadiusOfHNinteractions2((Xtotal-Xelastic)/pi/10.);
503 // SetProbabilityOfAnnihilation(1.);
504 // SetProbabilityOfAnnihilation(0.);
505 // //=======================================================
506 
507 //-----------------------------------------------------------------------------------
508 
509  SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
510  // (GeV/c)^(-2))
511 //G4cout<<"Slope "<<GetSlope()<<G4endl;
512 //-----------------------------------------------------------------------------------
513  SetGamma0( GetSlope()*Xtotal/10./2./pi );
514 
515 //----------- Parameters of elastic scattering --------------------------------------
516  // Gaussian parametrization of
517  // elastic scattering amplitude assumed
518  SetAvaragePt2ofElasticScattering(1./(Xtotal*Xtotal/16./pi/Xelastic/0.3894)*GeV*GeV);
519 //G4cout<<"AvaragePt2ofElasticScattering "<<GetAvaragePt2ofElasticScattering()<<G4endl;
520 //----------- Parameters of excitations ---------------------------------------------
521 
522  G4double Xinel=Xtotal-Xelastic; // Uzhi 25.04.2012
523 //G4cout<<"Param ProjectilePDGcode "<<ProjectilePDGcode<<G4endl;
524  if( ProjectilePDGcode > 1000 ) //------Projectile is baryon --------
525  {
526  SetMagQuarkExchange(1.84);//(3.63);
527  SetSlopeQuarkExchange(0.7);//(1.2);
528  SetDeltaProbAtQuarkExchange(0.);
529  if(NumberOfTargetNucleons > 26) {SetProbOfSameQuarkExchange(1.);}
530  else {SetProbOfSameQuarkExchange(0.);}
531 
532  SetProjMinDiffMass(1.16); // GeV
533  SetProjMinNonDiffMass(1.16); // GeV
534 // SetProbabilityOfProjDiff(0.805*std::exp(-0.35*Ylab)); // Uzhi 21.05.2012
535  SetProbabilityOfProjDiff(6./Xinel+1.5/ECMSsqr); // Uzhi 25.04.2012
536 
537  SetTarMinDiffMass(1.16); // GeV
538  SetTarMinNonDiffMass(1.16); // GeV
539 // SetProbabilityOfTarDiff(0.805*std::exp(-0.35*Ylab)); // Uzhi 21.05.2012
540  SetProbabilityOfTarDiff(6./Xinel+1.5/ECMSsqr); // Uzhi 25.04.2012
541 // SetAveragePt2(0.15); // 0.15 GeV^2
542  SetAveragePt2(0.3); // 0.30 GeV^2 Uzhi 21.05.2012
543 
544  SetProbLogDistr(0.5); // Uzhi 21.05.2012
545  }
546  else if( ProjectilePDGcode < -1000 ) //------Projectile is anti_baryon --------
547  {
548  SetMagQuarkExchange(0.);
549  SetSlopeQuarkExchange(0.);
550  SetDeltaProbAtQuarkExchange(0.);
551  SetProbOfSameQuarkExchange(0.);
552 
553  SetProjMinDiffMass(ProjectileMass+0.22); // GeV
554  SetProjMinNonDiffMass(ProjectileMass+0.22); // GeV
555 // SetProbabilityOfProjDiff(0.805*std::exp(-0.35*Ylab)); // Uzhi 21.05.2012
556  SetProbabilityOfProjDiff(6./Xinel+1.5/ECMSsqr); // Uzhi 25.04.2012
557 
558  SetTarMinDiffMass(TargetMass+0.22); // GeV
559  SetTarMinNonDiffMass(TargetMass+0.22); // GeV
560 // SetProbabilityOfTarDiff(0.805*std::exp(-0.35*Ylab)); // Uzhi 21.05.2012
561  SetProbabilityOfTarDiff(6./Xinel+1.5/ECMSsqr); // Uzhi 25.04.2012
562 
563  SetAveragePt2(0.3); // 0.15 GeV^2 // Uzhi 21.05.2012
564 
565  SetProbLogDistr(0.5); // Uzhi 21.05.2012
566  }
567  else if( ProjectileabsPDGcode == 211 ||
568  ProjectilePDGcode == 111) //------Projectile is Pion -----------
569  {
570  SetMagQuarkExchange(240.);
571  SetSlopeQuarkExchange(2.);
572  SetDeltaProbAtQuarkExchange(0.56); //(0.35);
573 
574  SetProjMinDiffMass(0.5); // GeV
575  SetProjMinNonDiffMass(0.5); // GeV 0.3
576 // SetProbabilityOfProjDiff(0.); // Uzhi 3.06.2012
577  SetProbabilityOfProjDiff((6.2-3.7*std::exp(-sqr(SqrtS-7.)/16.))/Xinel*0.);
578 
579  SetTarMinDiffMass(1.16); // GeV
580  SetTarMinNonDiffMass(1.16); // GeV
581 // SetProbabilityOfTarDiff(0.8*std::exp(-0.6*(Ylab-3.))); // Uzhi 3.06.2012
582  SetProbabilityOfTarDiff((2.+22./ECMSsqr)/Xinel);
583 
584  SetAveragePt2(0.3); // GeV^2 7 June 2011
585  SetProbLogDistr(0.); // Uzhi 21.05.2012
586 
587  SetProbLogDistr(1.);
588  }
589  else if( (ProjectileabsPDGcode == 321) ||
590  (ProjectileabsPDGcode == 311) ||
591  (ProjectilePDGcode == 130) ||
592  (ProjectilePDGcode == 310)) //Projectile is Kaon
593  {
594  SetMagQuarkExchange(40.);
595  SetSlopeQuarkExchange(2.25);
596  SetDeltaProbAtQuarkExchange(0.6);
597 
598  SetProjMinDiffMass(0.6); // GeV 0.7 0.6
599  SetProjMinNonDiffMass(0.6); // GeV 0.7 0.6
600 // SetProbabilityOfProjDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
601  SetProbabilityOfProjDiff(0.*4.7/Xinel); // Uzhi 5.06.2012
602 
603  SetTarMinDiffMass(1.1); // GeV
604  SetTarMinNonDiffMass(1.1); // GeV
605 // SetProbabilityOfTarDiff(0.45*std::pow(s/GeV/GeV,-0.5));// 40/32 X-dif/X-inel
606  SetProbabilityOfTarDiff(1.5/Xinel); // Uzhi 5.06.2012
607  SetAveragePt2(0.3); // GeV^2 7 June 2011
608  SetProbLogDistr(1.); // Uzhi 5.06.2012
609  }
610  else //------Projectile is undefined,
611  //------Nucleon assumed
612  {
613 /* // Uzhi 6.06.2012
614  SetMagQuarkExchange(1.85); // 7 June 2011
615  SetSlopeQuarkExchange(0.7); // 7 June 2011
616  SetDeltaProbAtQuarkExchange(0.); // 7 June 2011
617 
618  SetProjMinDiffMass((940.+160.*MeV)/GeV); // particle->GetPDGMass()
619  SetProjMinNonDiffMass((940.+160.*MeV)/GeV); // particle->GetPDGMass()
620  SetProbabilityOfProjDiff(0.805*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
621 
622  SetTarMinDiffMass(1.16); // GeV
623  SetTarMinNonDiffMass(1.16); // GeV
624  SetProbabilityOfTarDiff(0.805*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
625 */
626 
627  SetMagQuarkExchange(0.);
628  SetSlopeQuarkExchange(0.);
629  SetDeltaProbAtQuarkExchange(0.);
630  SetProbOfSameQuarkExchange(0.);
631 
632  SetProjMinDiffMass(ProjectileMass+0.22); // GeV
633  SetProjMinNonDiffMass(ProjectileMass+0.22); // GeV
634  SetProbabilityOfProjDiff(6./Xinel+1.5/ECMSsqr); // Uzhi 25.04.2012
635 
636  SetTarMinDiffMass(TargetMass+0.22); // GeV
637  SetTarMinNonDiffMass(TargetMass+0.22); // GeV
638  SetProbabilityOfTarDiff(6./Xinel+1.5/ECMSsqr); // Uzhi 25.04.2012
639 
640  SetAveragePt2(0.3); // 0.15 GeV^2 Uzhi 21.05.2012
641  SetProbLogDistr(0.5); // Uzhi 21.05.2012
642  }
643 
644 // if(theA > 4) SetProbabilityOfProjDiff(0.); // Uzhi 6.07.2012 Closed
645 
646 //G4cout<<"Param Get Min Dif "<<GetProjMinNonDiffMass()<<G4endl;
647 
648 // ---------- Set parameters of a string kink -------------------------------
649  SetPt2Kink(6.*GeV*GeV);
650  G4double Puubar(1./3.), Pddbar(1./3.), Pssbar(1./3.); // SU(3) symmetry
651 // G4double Puubar(0.41 ), Pddbar(0.41 ), Pssbar(0.18 ); // Broken SU(3) symmetry
652  SetQuarkProbabilitiesAtGluonSplitUp(Puubar, Pddbar, Pssbar);
653 
654 // --------- Set parameters of nuclear destruction--------------------
655 
656  if( ProjectileabsPDGcode < 1000 ) // Meson projectile
657  {
658  SetMaxNumberOfCollisions(Plab,2.); //3.); ##############################
659  SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/
660  (1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
661 
662  SetR2ofNuclearDestruction(1.5*fermi*fermi);
663 
664  SetDofNuclearDestruction(0.3);
665  SetPt2ofNuclearDestruction((0.035+0.04*std::exp(4.*(Ylab-2.5))/
666  (1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
667 //G4cout<<"Parm Pt2 Y "<<(0.035+0.04*std::exp(4.*(Ylab-2.5))/(1.+std::exp(4.*(Ylab-2.5))))<<" "<<Ylab<<G4endl;
668  SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
669 
670  SetExcitationEnergyPerWoundedNucleon(100.*MeV);
671  } else if( ProjectilePDGcode < -1000 ) // for anti-baryon projectile
672  {
673 //G4cout<<"Nucl destruct Anti Bar"<<G4endl;
674 
675  SetMaxNumberOfCollisions(Plab,2.); //3.); ##############################
676  SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/
677  (1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
678 
679  SetR2ofNuclearDestruction(1.5*fermi*fermi);
680 
681  SetDofNuclearDestruction(0.3);
682  SetPt2ofNuclearDestruction((0.035+0.04*std::exp(4.*(Ylab-2.5))/
683  (1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
684  SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
685 
686  SetExcitationEnergyPerWoundedNucleon(100.*MeV);
687 
688  if(Plab < 2.) // 2 GeV/c
689  { // For slow anti-baryon we have to garanty putting on mass-shell
690  SetCofNuclearDestruction(0.);
691  SetR2ofNuclearDestruction(1.5*fermi*fermi);
692  SetDofNuclearDestruction(0.01);
693  SetPt2ofNuclearDestruction(0.035*GeV*GeV);
694  SetMaxPt2ofNuclearDestruction(0.04*GeV*GeV);
695 // SetExcitationEnergyPerWoundedNucleon(0.); // ?????
696  }
697  } else // Projectile baryon assumed
698  {
699  SetMaxNumberOfCollisions(Plab,2.); //3.); ##############################
700  SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/
701  (1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
702 
703  SetR2ofNuclearDestruction(1.5*fermi*fermi);
704 
705  SetDofNuclearDestruction(0.3);
706  SetPt2ofNuclearDestruction((0.035+0.04*std::exp(4.*(Ylab-2.5))/
707  (1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
708  SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
709 
710  SetExcitationEnergyPerWoundedNucleon(100.*MeV);
711  }
712 
713 //SetCofNuclearDestruction(0.47*std::exp(2.*(Ylab-2.5))/(1.+std::exp(2.*(Ylab-2.5))));
714 //SetPt2ofNuclearDestruction((0.035+0.1*std::exp(4.*(Ylab-3.))/(1.+std::exp(4.*(Ylab-3.))))*GeV*GeV);
715 
716 //SetMagQuarkExchange(120.); // 210. PipP
717 //SetSlopeQuarkExchange(2.0);
718 //SetDeltaProbAtQuarkExchange(0.6);
719 //SetProjMinDiffMass(0.7); // GeV 1.1
720 //SetProjMinNonDiffMass(0.7); // GeV
721 //SetProbabilityOfProjDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
722 //SetTarMinDiffMass(1.1); // GeV
723 //SetTarMinNonDiffMass(1.1); // GeV
724 //SetProbabilityOfTarDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
725 //
726 //SetAveragePt2(0.3); // GeV^2
727 //------------------------------------
728 //SetProbabilityOfElasticScatt(1.,1.); //(Xtotal, Xelastic);
729 //SetProbabilityOfProjDiff(1.*0.62*std::pow(s/GeV/GeV,-0.51)); // 0->1
730 //SetProbabilityOfTarDiff(4.*0.62*std::pow(s/GeV/GeV,-0.51)); // 2->4
731 //SetAveragePt2(0.3); //(0.15);
732 //SetAvaragePt2ofElasticScattering(0.);
733 
734 //SetMaxNumberOfCollisions(Plab,6.); //(4.*(Plab+0.01),Plab); //6.); // ##########
735 //SetAveragePt2(0.15);
736 //G4cout<<"Cnd "<<GetCofNuclearDestruction()<<G4endl;
737 //SetCofNuclearDestruction(0.4);// (0.2); //(0.4);
738 //SetExcitationEnergyPerWoundedNucleon(0.*MeV); //(75.*MeV);
739 //SetDofNuclearDestruction(0.);
740 //SetPt2ofNuclearDestruction(0.*GeV*GeV); //(0.168*GeV*GeV);
741 //G4cout<<"Pt2 "<<GetPt2ofNuclearDestruction()/GeV/GeV<<G4endl;
742 //G4int Uzhi; G4cin>>Uzhi;
743 }
744 //**********************************************************************************************
Generated on Sat May 25 2013 14:34:05 for Geant4 by   doxygen 1.8.4