Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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
49G4FTFParameters::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
68G4FTFParameters::~G4FTFParameters()
69{}
70//**********************************************************************************************
71G4FTFParameters::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/*
290G4cout<<"Plab Xtotal, Xelastic Xinel Xftf "<<Plab<<" "<<Xtotal<<" "<<Xelastic<<" "<<Xtotal-Xelastic<<" "<<Xtotal-Xelastic-Xannihilation<<G4endl;
291G4cout<<"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/*
477G4cout<<"Plab Xtotal, Xelastic Xinel Xftf "<<Plab<<" "<<Xtotal<<" "<<Xelastic<<" "<<Xtotal-Xelastic<<" "<<Xtotal-Xelastic-Xannihilation<<G4endl;
478if(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//**********************************************************************************************
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4bool
bool G4bool
Definition: G4Types.hh:67
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4FTFParameters::SetTarMinNonDiffMass
void SetTarMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:273
G4FTFParameters::SetProjMinDiffMass
void SetProjMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:264
G4FTFParameters::SetTotalCrossSection
void SetTotalCrossSection(const G4double Xtotal)
Definition: G4FTFParameters.hh:218
G4FTFParameters::SetExcitationEnergyPerWoundedNucleon
void SetExcitationEnergyPerWoundedNucleon(const G4double aValue)
Definition: G4FTFParameters.hh:322
G4FTFParameters::FTFxsManager
G4ChipsComponentXS * FTFxsManager
Definition: G4FTFParameters.hh:164
G4FTFParameters::SetProbabilityOfElasticScatt
void SetProbabilityOfElasticScatt(const G4double Xtotal, const G4double Xelastic)
Definition: G4FTFParameters.hh:227
G4FTFParameters::SetElastisCrossSection
void SetElastisCrossSection(const G4double Xelastic)
Definition: G4FTFParameters.hh:221
G4FTFParameters::SetQuarkProbabilitiesAtGluonSplitUp
void SetQuarkProbabilitiesAtGluonSplitUp(const G4double Puubar, const G4double Pddbar, const G4double Pssbar)
Definition: G4FTFParameters.hh:288
G4FTFParameters::SetProbabilityOfTarDiff
void SetProbabilityOfTarDiff(const G4double aValue)
Definition: G4FTFParameters.hh:275
G4FTFParameters::SetMaxPt2ofNuclearDestruction
void SetMaxPt2ofNuclearDestruction(const G4double aValue)
Definition: G4FTFParameters.hh:329
G4FTFParameters::ProbOfSameQuarkExchange
G4double ProbOfSameQuarkExchange
Definition: G4FTFParameters.hh:181
G4FTFParameters::SetPt2Kink
void SetPt2Kink(const G4double aValue)
Definition: G4FTFParameters.hh:285
G4FTFParameters::SetSlope
void SetSlope(const G4double Slope)
Definition: G4FTFParameters.hh:243
G4FTFParameters::SetDeltaProbAtQuarkExchange
void SetDeltaProbAtQuarkExchange(const G4double aValue)
Definition: G4FTFParameters.hh:259
G4FTFParameters::SetMagQuarkExchange
void SetMagQuarkExchange(const G4double aValue)
Definition: G4FTFParameters.hh:255
G4FTFParameters::SetAvaragePt2ofElasticScattering
void SetAvaragePt2ofElasticScattering(const G4double aPt2)
Definition: G4FTFParameters.hh:250
G4FTFParameters::SetTarMinDiffMass
void SetTarMinDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:271
G4FTFParameters::SetProjMinNonDiffMass
void SetProjMinNonDiffMass(const G4double aValue)
Definition: G4FTFParameters.hh:266
G4FTFParameters::SetR2ofNuclearDestruction
void SetR2ofNuclearDestruction(const G4double aValue)
Definition: G4FTFParameters.hh:319
G4FTFParameters::SetGamma0
void SetGamma0(const G4double Gamma0)
Definition: G4FTFParameters.hh:246
G4FTFParameters::SetAveragePt2
void SetAveragePt2(const G4double aValue)
Definition: G4FTFParameters.hh:278
G4FTFParameters::SetProbabilityOfAnnihilation
void SetProbabilityOfAnnihilation(const G4double aValue)
Definition: G4FTFParameters.hh:237
G4FTFParameters::SetDofNuclearDestruction
void SetDofNuclearDestruction(const G4double aValue)
Definition: G4FTFParameters.hh:325
G4FTFParameters::SetRadiusOfHNinteractions2
void SetRadiusOfHNinteractions2(const G4double Radius2)
Definition: G4FTFParameters.hh:240
G4FTFParameters::SetPt2ofNuclearDestruction
void SetPt2ofNuclearDestruction(const G4double aValue)
Definition: G4FTFParameters.hh:327
G4FTFParameters::SetProbLogDistr
void SetProbLogDistr(const G4double aValue)
Definition: G4FTFParameters.hh:281
G4FTFParameters::SetCofNuclearDestruction
void SetCofNuclearDestruction(const G4double aValue)
Definition: G4FTFParameters.hh:317
G4FTFParameters::SetMaxNumberOfCollisions
void SetMaxNumberOfCollisions(const G4double aValue, const G4double bValue)
Definition: G4FTFParameters.hh:299
G4FTFParameters::SetSlopeQuarkExchange
void SetSlopeQuarkExchange(const G4double aValue)
Definition: G4FTFParameters.hh:257
G4FTFParameters::SetProbOfSameQuarkExchange
void SetProbOfSameQuarkExchange(const G4double aValue)
Definition: G4FTFParameters.hh:261
G4FTFParameters::FTFXannihilation
G4double FTFXannihilation
Definition: G4FTFParameters.hh:169
G4FTFParameters::SetProbabilityOfProjDiff
void SetProbabilityOfProjDiff(const G4double aValue)
Definition: G4FTFParameters.hh:268
G4FTFParameters::SetInelasticCrossSection
void SetInelasticCrossSection(const G4double Xinelastic)
Definition: G4FTFParameters.hh:224
G4KaonMinus::KaonMinus
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:119
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
sqr
T sqr(const T &x)
Definition: templates.hh:145