| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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
66 // Note that the very first tune (with index
67 // set of parameters, which does not need to
68 // the for loop below starts from 1 and not
69 // The check whether an alternative tune has
70 // level of the G4FTFParamCollection::SetTun
71 for ( G4int indexTune = 1; indexTune < G4FTF
72 fArrayParCollBaryonProj[indexTune].SetTune
73 fArrayParCollMesonProj[indexTune].SetTune(
74 fArrayParCollPionProj[indexTune].SetTune(i
75 }
76
77 StringMass = new G4LundStringFragmentation;
78 Reset();
79 csGGinstance =
80 G4CrossSectionDataSetRegistry::Instance()-
81 if (!csGGinstance) {
82 csGGinstance = new G4ComponentGGHadronNucl
83 }
84
85 EnableDiffDissociationForBGreater10 = G4Hadr
86
87 // Set parameters of a string kink
88 SetPt2Kink( 0.0*GeV*GeV ); // To switch off
89 G4double Puubar( 1.0/3.0 ), Pddbar( 1.0/3.0
90 //G4double Puubar( 0.41 ), Pddbar( 0.41 ), P
91 SetQuarkProbabilitiesAtGluonSplitUp( Puubar,
92 }
93
94 //============================================
95
96 void G4FTFParameters::InitForInteraction( cons
97 G4in
98 {
99 Reset();
100
101 G4int ProjectilePDGcode = particle->Ge
102 G4int ProjectileabsPDGcode = std::abs( Pr
103 G4double ProjectileMass = particle->Ge
104 G4double ProjectileMass2 = ProjectileMa
105
106 G4int ProjectileBaryonNumber( 0 ), AbsProjec
107 G4bool ProjectileIsNucleus = false;
108
109 if ( std::abs( particle->GetBaryonNumber() )
110 ProjectileIsNucleus = true;
111 ProjectileBaryonNumber = particle->GetB
112 AbsProjectileBaryonNumber = std::abs( Proj
113 AbsProjectileCharge = std::abs( G4in
114 if ( ProjectileBaryonNumber > 1 ) {
115 ProjectilePDGcode = 2212; ProjectileabsP
116 } else {
117 ProjectilePDGcode = -2212; Projectileabs
118 }
119 ProjectileMass = G4Proton::Proton()->GetP
120 ProjectileMass2 = sqr( ProjectileMass );
121 }
122
123 G4double TargetMass = G4Proton::Proton()->G
124 G4double TargetMass2 = TargetMass * TargetMa
125
126 G4double Plab = PlabPerParticle;
127 G4double Elab = std::sqrt( Plab*Plab + Proje
128 G4double KineticEnergy = Elab - ProjectileMa
129
130 G4double S = ProjectileMass2 + TargetMass2 +
131
132 #ifdef debugFTFparams
133 G4cout << "--------- FTF Parameters --------
134 << ProjectilePDGcode << " " << Plab <
135 << " " << KineticEnergy << G4endl <<
136 #endif
137
138 G4double Ylab, Xtotal( 0.0 ), Xelastic( 0.0
139 G4int NumberOfTargetNucleons;
140
141 Ylab = 0.5 * G4Log( (Elab + Plab)/(Elab - Pl
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 /= (G
151 ProjectileMass /= GeV; ProjectileMass2 /= (G
152
153 Plab /= GeV;
154 G4double Xftf = 0.0;
155
156 G4int NumberOfTargetProtons = theZ;
157 G4int NumberOfTargetNeutrons = theA - theZ;
158 NumberOfTargetNucleons = NumberOfTargetProto
159
160 // ---------- hadron projectile ------------
161 if ( AbsProjectileBaryonNumber <= 1 ) { //
162
163 // Interaction on P
164 G4double xTtP = csGGinstance->GetTotalIsot
165 G4double xElP = csGGinstance->GetElasticIs
166
167 // Interaction on N
168 G4double xTtN = csGGinstance->GetTotalIsot
169 G4double xElN = csGGinstance->GetElasticIs
170
171 // Average properties of h+N interactions
172 Xtotal = ( NumberOfTargetProtons * xTtP
173 Xelastic = ( NumberOfTargetProtons * xElP
174 Xannihilation = 0.0;
175
176 Xtotal /= millibarn;
177 Xelastic /= millibarn;
178
179 #ifdef debugFTFparams
180 G4cout<<"Estimated cross sections (total a
181 #endif
182 }
183
184 // ---------- nucleus projectile -----------
185 if ( ProjectileIsNucleus && ProjectileBary
186
187 #ifdef debugFTFparams
188 G4cout<<"Projectile is a nucleus: A and Z
189 #endif
190
191 const G4ParticleDefinition* Proton = G4Pro
192 // Interaction on P
193 G4double XtotPP = csGGinstance->GetTotalIs
194 G4double XelPP = csGGinstance->GetElastic
195
196 const G4ParticleDefinition* Neutron = G4Ne
197 // Interaction on N
198 G4double XtotPN = csGGinstance->GetTotalIs
199 G4double XelPN = csGGinstance->GetElastic
200
201 #ifdef debugFTFparams
202 G4cout << "XsPP (total and elastic) " << X
203 << "XsPN (total and elastic) " << X
204 #endif
205
206 Xtotal = (
207 AbsProjectileCharge * Number
208 ( AbsProjectileBaryonNumber
209 NumberOfTargetNeutrons *
210 +
211 ( AbsProjectileCharge * Numb
212 ( AbsProjectileBaryonNumbe
213 NumberOfTargetProtons
214 ) / ( AbsProjectileBaryonNumbe
215 Xelastic= (
216 AbsProjectileCharge * Number
217 ( AbsProjectileBaryonNumber
218 NumberOfTargetNeutrons *
219 +
220 ( AbsProjectileCharge * Numb
221 ( AbsProjectileBaryonNumbe
222 NumberOfTargetProtons
223 ) / ( AbsProjectileBaryonNumbe
224
225 Xannihilation = 0.0;
226 Xtotal /= millibarn;
227 Xelastic /= millibarn;
228 }
229
230 // ---------- The projectile is anti-baryon
231 // anti Sigma^0_c
232 if ( ProjectilePDGcode >= -4112 && Project
233 // Only non-strange and strange baryons ar
234
235 #ifdef debugFTFparams
236 G4cout<<"Projectile is a anti-baryon or an
237 G4cout<<"(Only non-strange and strange bar
238 #endif
239
240 G4double X_a( 0.0 ), X_b( 0.0 ), X_c( 0.0
241 G4double MesonProdThreshold = ProjectileMa
242 ( 2.0 * 0.14
243
244 if ( PlabPerParticle < 40.0*MeV ) { // Low
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 sectio
252 G4double LogS = G4Log( ECMSsqr / 33.0625
253 G4double Xasmpt = 36.04 + 0.304*LogS*Log
254 LogS = G4Log( SqrtS / 20.74 );
255 G4double Basmpt = 11.92 + 0.3036*LogS*Lo
256 G4double R0 = std::sqrt( 0.40874044*Xasm
257
258 G4double FlowF = SqrtS / std::sqrt( ECMS
259 Targ
260 - 2.
261 - 2.
262
263 Xtotal = Xasmpt * ( 1.0 + 13.55*FlowF/R0
264 (1.0 - 4.47/Sq
265
266 Xasmpt = 4.4 + 0.101*LogS*LogS; // mb
267 Xelastic = Xasmpt * ( 1.0 + 59.27*FlowF/
268 (1.0 - 6.95/
269
270 X_a = 25.0*FlowF; // mb, 3-shirts diagr
271
272 if ( SqrtS < MesonProdThreshold ) {
273 X_b = 3.13 + 140.0*G4Pow::GetInstance(
274 Xelastic -= 3.0*X_b; // Xel-X(PbarP->
275 } else {
276 X_b = 6.8/SqrtS; // mb anti-quark-qua
277 Xelastic -= 3.0*X_b; // Xel-X(PbarP->
278 }
279
280 X_c = 2.0*FlowF*sqr( ProjectileMass + Ta
281
282 X_d = 23.3/ECMSsqr; // mb anti-quark-qu
283 }
284
285 G4double Xann_on_P( 0.0), Xann_on_N( 0.0 )
286
287 if ( ProjectilePDGcode == -2212 ) {
288 Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_
289 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_
290 } else if ( ProjectilePDGcode == -2112 ) {
291 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_
292 Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_
293 } else if ( ProjectilePDGcode == -3122 ) {
294 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_
295 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_
296 } else if ( ProjectilePDGcode == -3112 ) {
297 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_
298 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_
299 } else if ( ProjectilePDGcode == -3212 ) {
300 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_
301 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_
302 } else if ( ProjectilePDGcode == -3222 ) {
303 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_
304 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_
305 } else if ( ProjectilePDGcode == -3312 ) {
306 Xann_on_P = X_a + X_b*1.0 + X_c*1.0 + X_
307 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_
308 } else if ( ProjectilePDGcode == -3322 ) {
309 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_
310 Xann_on_N = X_a + X_b*1.0 + X_c*1.0 + X_
311 } else if ( ProjectilePDGcode == -3334 ) {
312 Xann_on_P = X_a + X_b*0.0 + X_c*0.0 + X_
313 Xann_on_N = X_a + X_b*0.0 + X_c*0.0 + X_
314 } else {
315 G4cout << "Unknown anti-baryon for FTF a
316 }
317
318 //G4cout << "Sum " << Xann_on_P <
319
320 if ( ! ProjectileIsNucleus ) { // Project
321 Xannihilation = ( NumberOfTargetProtons
322 / NumberOfTargetNucleons
323 } else { // Projectile is a nucleus
324 Xannihilation = (
325 ( AbsProjectileCharge
326 ( AbsProjectileBaryo
327 NumberOfTargetNeutro
328 +
329 ( AbsProjectileCharge
330 ( AbsProjectileBaryo
331 NumberOfTargetProton
332 ) / ( AbsProjectileBaryo
333 }
334
335 //G4double Xftf = 0.0;
336 MesonProdThreshold = ProjectileMass + Targ
337 if ( SqrtS > MesonProdThreshold ) {
338 Xftf = 36.0 * ( 1.0 - MesonProdThreshold
339 }
340
341 Xtotal = Xelastic + Xannihilation + Xftf;
342
343 #ifdef debugFTFparams
344 G4cout << "Plab Xtotal, Xelastic Xinel Xf
345 << " " << Xtotal - Xelastic << " "
346 << "Plab Xelastic/Xtotal, Xann/Xin
347 << Xannihilation/(Xtotal - Xelastic
348 #endif
349
350 }
351
352 if ( Xtotal == 0.0 ) { // Projectile is und
353
354 const G4ParticleDefinition* Proton = G4Pro
355 // Interaction on P
356 G4double XtotPP = csGGinstance->GetTotalIs
357 G4double XelPP = csGGinstance->GetElastic
358
359 // Interaction on N
360 G4double XtotPN = csGGinstance->GetTotalIs
361 G4double XelPN = csGGinstance->GetElastic
362
363 Xtotal = ( NumberOfTargetProtons * Xtot
364 / NumberOfTargetNucleons;
365 Xelastic = ( NumberOfTargetProtons * XelP
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 c
378 SetProbabilityOfElasticScatt( Xtotal, Xelast
379
380 SetRadiusOfHNinteractions2( Xtotal/pi/10.0 )
381
382 if ( ( Xtotal - Xelastic ) == 0.0 ) {
383 SetProbabilityOfAnnihilation( 0.0 );
384 } else {
385 SetProbabilityOfAnnihilation( Xannihilatio
386 }
387
388 if(Xelastic > 0.0) {
389 SetSlope( Xtotal*Xtotal/16.0/pi/Xelastic/0
390
391 // Parameters of elastic scattering
392 // Gaussian parametrization of elastic sca
393 SetAvaragePt2ofElasticScattering( 1.0/( Xt
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" <<
404 G4cout << "AvaragePt2ofElasticScattering " <
405 G4cout<<"Parameters of excitation for projec
406 #endif
407
408 if ( (ProjectilePDGcode == 2212) || (Project
409
410 const G4int indexTune = G4FTFTunings::Inst
411
412 // A process probability is parameterized
413 // y is a rapidity of a partcle in the tar
414
415 // Proc# A1 B1 A2
416 /* original hadr-string-diff-V10-03-07 (si
417 SetParams( 0, 13.71, 1.75, -2
418 SetParams( 1, 25.0, 1.0, -5
419 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*1
420 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*1
421 SetParams( 4, 1.0, 0.0 , -2
422 */
423
424 // Proc#
425 SetParams( 0, fArrayParCollBaryonProj[inde
426 fArrayParCollBaryonProj[indexTune]
427 fArrayParCollBaryonProj[indexTune].GetPr
428 fArrayParCollBaryonProj[indexTune]
429 fArrayParCollBaryonProj[indexTune].GetPr
430 fArrayParCollBaryonProj[indexTune]
431 fArrayParCollBaryonProj[indexTune].GetPr
432 SetParams( 1, fArrayParCollBaryonProj[inde
433 fArrayParCollBaryonProj[indexTune]
434 fArrayParCollBaryonProj[indexTune].GetPr
435 fArrayParCollBaryonProj[indexTune]
436 fArrayParCollBaryonProj[indexTune].GetPr
437 fArrayParCollBaryonProj[indexTune]
438 fArrayParCollBaryonProj[indexTune].GetPr
439 if ( Xinel > 0.0 ) {
440 SetParams( 2, 6.0/Xinel, 0.0, -6.0/Xinel
441 SetParams( 3, 6.0/Xinel, 0.0, -6.0/Xinel
442
443 SetParams( 4, fArrayParCollBaryonProj[in
444 fArrayParCollBaryonProj[indexTune].Get
445 fArrayParCollBaryonProj[indexTune].Get
446 fArrayParCollBaryonProj[indexTune].Get
447 fArrayParCollBaryonProj[indexTune].Get
448 fArrayParCollBaryonProj[indexTune].Get
449 fArrayParCollBaryonProj[indexTune].Get
450 } else { // if Xinel=0., zero everything
451 SetParams( 2, 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
453 SetParams( 4, 0.0, 0.0, 0.0, 0.0, 0.0, 0
454 }
455
456 if ( (AbsProjectileBaryonNumber > 10 || Nu
457 // It is not decided what to do with dif
458 // For the moment both ProjDiffDisso & T
459 // so both projectile and target diffrac
460 if ( ! fArrayParCollBaryonProj[indexTune
461 SetParams( 2, 0.0, 0.0, 0.0, 0.0, 0.0
462 if ( ! fArrayParCollBaryonProj[indexTune
463 SetParams( 3, 0.0, 0.0, 0.0, 0.0, 0.0
464 }
465
466 SetDeltaProbAtQuarkExchange( fArrayParColl
467
468 if ( NumberOfTargetNucleons > 26 ) {
469 SetProbOfSameQuarkExchange( 1.0 );
470 } else {
471 SetProbOfSameQuarkExchange( fArrayParCol
472 }
473
474 SetProjMinDiffMass( fArrayParCollBaryon
475 SetProjMinNonDiffMass( fArrayParCollBaryon
476
477 SetTarMinDiffMass( fArrayParCollBaryon
478 SetTarMinNonDiffMass( fArrayParCollBaryon
479
480 SetAveragePt2( fArrayParCollBaryon
481 SetProbLogDistrPrD( fArrayParCollBaryon
482 SetProbLogDistr( fArrayParCollBaryon
483
484 } else if ( ProjectilePDGcode == -2212 ||
485
486 // Below, in the call to the G4FTFTunings:
487 // as projectile, instead of the real one,
488 // we assume the same treatment for anti_p
489 // whereas all other parameters for anti_p
490 const G4int indexTune = G4FTFTunings::Inst
491
492 // Proc# A1 B1 A2
493 SetParams( 0, 0.0 , 0.0 , 0
494 SetParams( 1, 0.0 , 0.0 , 0
495 if ( Xinel > 0.) {
496 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel
497 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel
498 SetParams( 4, 1.0, 0.0 ,
499 } else {
500 SetParams( 2, 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
502 SetParams( 4, 0.0, 0.0, 0.0, 0.0, 0.0, 0
503 }
504
505 if ( AbsProjectileBaryonNumber > 10 || N
506 // It is not decided what to do with dif
507 // For the moment both ProjDiffDisso & T
508 // so both projectile and target diffrac
509 if ( ! fArrayParCollBaryonProj[indexTune
510 SetParams( 2, 0.0, 0.0 ,
511 if ( ! fArrayParCollBaryonProj[indexTune
512 SetParams( 3, 0.0, 0.0 ,
513 }
514
515 SetDeltaProbAtQuarkExchange( 0.0 );
516 SetProbOfSameQuarkExchange( 0.0 );
517 SetProjMinDiffMass( ProjectileMass + 0.22
518 SetProjMinNonDiffMass( ProjectileMass + 0.
519 SetTarMinDiffMass( TargetMass + 0.22 );
520 SetTarMinNonDiffMass( TargetMass + 0.22 );
521 SetAveragePt2( 0.3 );
522 SetProbLogDistrPrD( 0.55 );
523 SetProbLogDistr( 0.55 );
524
525 } else if ( ProjectileabsPDGcode == 211 ||
526
527 const G4int indexTune = G4FTFTunings::Inst
528
529 // Proc# A1 B1 A2
530 /* --> original code
531 SetParams( 0, 150.0, 1.8 , -247.
532 SetParams( 1, 5.77, 0.6 , -5.7
533 SetParams( 2, 2.27, 0.5 , -98052.
534 SetParams( 3, 7.0, 0.9, -85.2
535 SetParams( 4, 1.0, 0.0 , -11.0
536 */
537 // Proc#
538 SetParams( 0, fArrayParCollPionProj[indexT
539 fArrayParCollPionProj[indexTune].G
540 fArrayParCollPionProj[indexTune].GetProc
541 fArrayParCollPionProj[indexTune].G
542 fArrayParCollPionProj[indexTune].GetProc
543 fArrayParCollPionProj[indexTune].G
544 fArrayParCollPionProj[indexTune].GetProc
545 SetParams( 1, fArrayParCollPionProj[indexT
546 fArrayParCollPionProj[indexTune].G
547 fArrayParCollPionProj[indexTune].GetProc
548 fArrayParCollPionProj[indexTune].G
549 fArrayParCollPionProj[indexTune].GetProc
550 fArrayParCollPionProj[indexTune].G
551 fArrayParCollPionProj[indexTune].GetProc
552 SetParams( 2, fArrayParCollPionProj[indexT
553 fArrayParCollPionProj[indexTune].G
554 fArrayParCollPionProj[indexTune].GetProc
555 fArrayParCollPionProj[indexTune].G
556 fArrayParCollPionProj[indexTune].GetProc
557 fArrayParCollPionProj[indexTune].G
558 fArrayParCollPionProj[indexTune].GetProc
559 SetParams( 3, fArrayParCollPionProj[indexT
560 fArrayParCollPionProj[indexTune].G
561 fArrayParCollPionProj[indexTune].GetProc
562 fArrayParCollPionProj[indexTune].G
563 fArrayParCollPionProj[indexTune].GetProc
564 fArrayParCollPionProj[indexTune].G
565 fArrayParCollPionProj[indexTune].GetProc
566 SetParams( 4, fArrayParCollPionProj[indexT
567 fArrayParCollPionProj[indexTune].G
568 fArrayParCollPionProj[indexTune].GetProc
569 fArrayParCollPionProj[indexTune].G
570 fArrayParCollPionProj[indexTune].GetProc
571 fArrayParCollPionProj[indexTune].G
572 fArrayParCollPionProj[indexTune].GetProc
573
574 // NOTE: how can it be |ProjectileBaryonNu
575 //
576 if ( AbsProjectileBaryonNumber > 10 || N
577 if ( ! fArrayParCollPionProj[indexTune]
578 SetParams( 2, 0.0 , 0.0 ,
579 if ( ! fArrayParCollPionProj[indexTune]
580 SetParams( 3, 0.0 , 0.0 ,
581 }
582
583 /* original code -->
584 SetDeltaProbAtQuarkExchange( 0.56 );
585 SetProjMinDiffMass( 1.0 ); // G
586 SetProjMinNonDiffMass( 1.0 ); // G
587 SetTarMinDiffMass( 1.16 ); // G
588 SetTarMinNonDiffMass( 1.16 ); // G
589 SetAveragePt2( 0.3 ); // G
590 SetProbLogDistrPrD( 0.55 );
591 SetProbLogDistr( 0.55 );
592 */
593
594 // JVY update, Aug.8, 2018 --> Feb.14, 201
595 //
596 SetDeltaProbAtQuarkExchange( fArrayParColl
597 SetProjMinDiffMass( fArrayParCollPionProj[
598 SetProjMinNonDiffMass( fArrayParCollPionPr
599 SetTarMinDiffMass( fArrayParCollPionProj[i
600 SetTarMinNonDiffMass( fArrayParCollPionPro
601 SetAveragePt2( fArrayParCollPionProj[index
602 SetProbLogDistrPrD( fArrayParCollPionProj[
603 SetProbLogDistr( fArrayParCollPionProj[ind
604
605 // ---> end update
606
607 } else if ( ProjectileabsPDGcode == 321 ||
608 ProjectilePDGcode == 130 ||
609
610 // Proc# A1 B1 A2
611 SetParams( 0, 60.0 , 2.5 , 0
612 SetParams( 1, 6.0 , 1.0 , -24
613 SetParams( 2, 2.76, 1.2 , -22
614 SetParams( 3, 1.09, 0.5 , -8
615 SetParams( 4, 1.0, 0.0 , 0
616 if ( AbsProjectileBaryonNumber > 10 || N
617 SetParams( 2, 0.0 , 0.0 ,
618 SetParams( 3, 0.0 , 0.0 ,
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
631
632 if ( ProjectileabsPDGcode > 1000 ) { // T
633 // Proc# A1 B1
634 SetParams( 0, 13.71, 1.75,
635 SetParams( 1, 25.0, 1.0, -
636 if ( Xinel > 0.) {
637 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xin
638 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xin
639 SetParams( 4, 1.0, 0.0 ,
640 } else {
641 SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0
642 SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0
643 SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0
644 }
645
646 } else { // The projectile is a meson as
647 // Proc# A1 B1
648 SetParams( 0, 60.0 , 2.5 ,
649 SetParams( 1, 6.0 , 1.0 , -
650 SetParams( 2, 2.76, 1.2 , -
651 SetParams( 3, 1.09, 0.5 ,
652 SetParams( 4, 1.0, 0.0 ,
653 }
654
655 if ( AbsProjectileBaryonNumber > 10 || N
656 SetParams( 2, 0.0 , 0.0 ,
657 SetParams( 3, 0.0 , 0.0 ,
658 }
659
660 SetDeltaProbAtQuarkExchange( 0.0 );
661 SetProbOfSameQuarkExchange( 0.0 );
662
663 SetProjMinDiffMass( GetMinMass(particle
664 SetProjMinNonDiffMass( GetMinMass(particle
665
666 const G4ParticleDefinition* Neutron = G4Ne
667 SetTarMinDiffMass( GetMinMass(Neutron)/
668 SetTarMinNonDiffMass( GetMinMass(Neutron)/
669
670 SetAveragePt2( 0.3 ); // GeV^2
671 SetProbLogDistrPrD( 0.55 );
672 SetProbLogDistr( 0.55 );
673
674 }
675
676 #ifdef debugFTFparams
677 G4cout<<"DeltaProbAtQuarkExchange "<< GetDel
678 G4cout<<"ProbOfSameQuarkExchange "<< GetPro
679 G4cout<<"ProjMinDiffMass "<< GetPro
680 G4cout<<"ProjMinNonDiffMass "<< GetPro
681 G4cout<<"TarMinDiffMass "<< GetTar
682 G4cout<<"TarMinNonDiffMass "<< GetTar
683 G4cout<<"AveragePt2 "<< GetAve
684 G4cout<<"ProbLogDistrPrD "<< GetPro
685 G4cout<<"ProbLogDistrTrD "<< GetPro
686 #endif
687
688 // Set parameters of nuclear destruction
689
690 if ( ProjectileabsPDGcode < 1000 ) { // Mes
691
692 const G4int indexTune = G4FTFTunings::Inst
693
694 SetMaxNumberOfCollisions( Plab, 2.0 ); //
695 //
696 // target destruction
697 //
698 /* original code --->
699 SetCofNuclearDestruction( 0.00481*G4double
700 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 +
701
702 SetR2ofNuclearDestruction( 1.5*fermi*fermi
703 SetDofNuclearDestruction( 0.3 );
704 SetPt2ofNuclearDestruction( ( 0.035 + 0.04
705 ( 1.0
706 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV
707 SetExcitationEnergyPerWoundedNucleon( 40.0
708 */
709 double coeff = fArrayParCollMesonProj[inde
710 //
711 // NOTE (JVY): Set this switch to false/tr
712 //
713 if ( fArrayParCollMesonProj[indexTune].IsN
714 {
715 coeff *= G4double(NumberOfTargetNucleons
716 }
717 double exfactor = G4Exp( fArrayParCollMeso
718 * (Ylab-fArrayParCo
719 coeff *= exfactor;
720 coeff /= ( 1.+ exfactor );
721
722 SetCofNuclearDestruction( coeff );
723
724 SetR2ofNuclearDestruction( fArrayParCollMe
725 SetDofNuclearDestruction( fArrayParCollMes
726 coeff = fArrayParCollMesonProj[indexTune].
727 exfactor = G4Exp( fArrayParCollMesonProj[i
728 * (Ylab-fArrayParCollMeson
729 coeff *= exfactor;
730 coeff /= ( 1. + exfactor );
731 SetPt2ofNuclearDestruction( (fArrayParColl
732
733 SetMaxPt2ofNuclearDestruction( fArrayParCo
734 SetExcitationEnergyPerWoundedNucleon( fArr
735
736 } else if ( ProjectilePDGcode == -2212 ||
737
738 SetMaxNumberOfCollisions( Plab, 2.0 );
739
740 SetCofNuclearDestruction( 0.00481*G4double
741 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G
742 SetR2ofNuclearDestruction( 1.5*fermi*fermi
743 SetDofNuclearDestruction( 0.3 );
744 SetPt2ofNuclearDestruction( ( 0.035 + 0.04
745 ( 1.0
746 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV
747 SetExcitationEnergyPerWoundedNucleon( 40.0
748 if ( Plab < 2.0 ) { // 2 GeV/c
749 // For slow anti-baryon we have to garan
750 SetCofNuclearDestruction( 0.0 );
751 SetR2ofNuclearDestruction( 1.5*fermi*fer
752
753 SetDofNuclearDestruction( 0.01 );
754 SetPt2ofNuclearDestruction( 0.035*GeV*Ge
755 SetMaxPt2ofNuclearDestruction( 0.04*GeV*
756 }
757
758 } else { // Projectile baryon assumed
759
760 // Below, in the call to the G4FTFTunings:
761 // as projectile, instead of the real one,
762 // destruction, we assume for this categor
763 // as for "baryon".
764 const G4int indexTune = G4FTFTunings::Inst
765
766 // NOTE (JVY) FIXME !!! Will decide later
767 //
768 SetMaxNumberOfCollisions( Plab, 2.0 );
769
770 // projectile destruction - does NOT reall
771 //
772 double coeff = 0.;
773 coeff = fArrayParCollBaryonProj[indexTune]
774 //
775 // NOTE (JVY): Set this switch to false/tr
776 //
777 if ( fArrayParCollBaryonProj[indexTune].Is
778 {
779 coeff *= G4double(AbsProjectileBaryonNum
780 }
781 double exfactor = G4Exp( fArrayParCollBary
782 (Ylab-fArrayParCollBaryonProj[index
783 coeff *= exfactor;
784 coeff /= ( 1.+ exfactor );
785 SetCofNuclearDestructionPr( coeff );
786
787 // target desctruction
788 //
789 coeff = fArrayParCollBaryonProj[indexTune]
790 //
791 // NOTE (JVY): Set this switch to false/tr
792 //
793 if ( fArrayParCollBaryonProj[indexTune].Is
794 {
795 coeff *= G4double(NumberOfTargetNucleons
796 }
797 exfactor = G4Exp( fArrayParCollBaryonProj[
798 (Ylab-fArrayParCollBaryonProj[indexT
799 coeff *= exfactor;
800 coeff /= ( 1.+ exfactor );
801 SetCofNuclearDestruction( coeff );
802
803 SetR2ofNuclearDestruction( fArrayParCollBa
804 SetDofNuclearDestruction( fArrayParCollBar
805
806 coeff = fArrayParCollBaryonProj[indexTune]
807 exfactor = G4Exp( fArrayParCollBaryonProj[
808 (Ylab-fArrayParCollBaryonProj[indexT
809 coeff *= exfactor;
810 coeff /= ( 1. + exfactor );
811 SetPt2ofNuclearDestruction( (fArrayParColl
812
813 SetMaxPt2ofNuclearDestruction( fArrayParCo
814 SetExcitationEnergyPerWoundedNucleon( fArr
815
816 }
817
818 #ifdef debugFTFparams
819 G4cout<<"CofNuclearDestructionPr "
820 G4cout<<"CofNuclearDestructionTr "
821 G4cout<<"R2ofNuclearDestruction "
822 G4cout<<"DofNuclearDestruction "
823 G4cout<<"Pt2ofNuclearDestruction "
824 G4cout<<"ExcitationEnergyPerWoundedNucleon "
825 #endif
826
827 //SetCofNuclearDestruction( 0.47*G4Exp( 2.0*
828 //SetPt2ofNuclearDestruction( ( 0.035 + 0.1*
829
830 //SetMagQuarkExchange( 120.0 ); // 210
831 //SetSlopeQuarkExchange( 2.0 );
832 //SetDeltaProbAtQuarkExchange( 0.6 );
833 //SetProjMinDiffMass( 0.7 ); // GeV
834 //SetProjMinNonDiffMass( 0.7 ); // GeV
835 //SetProbabilityOfProjDiff( 0.0); // 0.8
836 //SetTarMinDiffMass( 1.1 ); // GeV
837 //SetTarMinNonDiffMass( 1.1 ); // GeV
838 //SetProbabilityOfTarDiff( 0.0 ); // 0.8
839
840 //SetAveragePt2( 0.0 ); // GeV
841 //------------------------------------
842 //SetProbabilityOfElasticScatt( 1.0, 1.0);
843 //SetProbabilityOfProjDiff( 1.0*0.62*G4Pow::
844 //SetProbabilityOfTarDiff( 4.0*0.62*G4Pow::G
845 //SetAveragePt2( 0.3 );
846 //SetAvaragePt2ofElasticScattering( 0.0 );
847
848 //SetMaxNumberOfCollisions( Plab, 6.0 ); //(
849 //SetAveragePt2( 0.15 );
850 //SetCofNuclearDestruction(-1.);//( 0.75 );
851 //SetExcitationEnergyPerWoundedNucleon(0.);/
852 //SetDofNuclearDestruction(0.);//( 0.2 ); //
853
854 /*
855 SetAveragePt2(0.3);
856 SetCofNuclearDestructionPr(0.);
857 SetCofNuclearDestruction(0.); //
858 SetExcitationEnergyPerWoundedNucleon(0.); //
859 SetDofNuclearDestruction(0.); //
860 SetPt2ofNuclearDestruction(0.); //
861 */
862
863 //SetExcitationEnergyPerWoundedNucleon(0.001
864 //SetPt2Kink( 0.0*GeV*GeV );
865
866 //SetRadiusOfHNinteractions2( Xtotal/pi/10.0
867 //SetRadiusOfHNinteractions2( (Xtotal - Xela
868 //SetProbabilityOfElasticScatt( 1.0, 0.0);
869 /*
870 G4cout << "Pt2 " << GetAveragePt2()<<" "<<Ge
871 G4cout << "Cnd " << GetCofNuclearDestruction
872 G4cout << "Dnd " << GetDofNuclearDestruction
873 G4cout << "Pt2 " << GetPt2ofNuclearDestructi
874 */
875
876 }
877
878 //============================================
879
880 G4double G4FTFParameters::GetMinMass( const G4
881 // The code is used for estimating the mini
882 // The indices used for minMassQDiQStr must
883 // quarks: d, u, s, c and b; enforcing this
884 G4double EstimatedMass = 0.0;
885 G4int partID = std::abs(aParticle->GetPDGEn
886 G4int Qleft = std::max( partID/100, 1
887 G4int Qright = std::max( (partID/ 10)%10, 1
888 if ( Qleft < 6 && Qright < 6 ) {
889 EstimatedMass = StringMass->minMassQQbarS
890 } else if ( Qleft < 6 && Qright > 6 ) {
891 G4int q1 = std::max( std::min( Qright/10,
892 G4int q2 = std::max( std::min( Qright%10,
893 EstimatedMass = StringMass->minMassQDiQSt
894 } else if ( Qleft > 6 && Qright < 6 ) {
895 G4int q1 = std::max( std::min( Qleft/10,
896 G4int q2 = std::max( std::min( Qleft%10,
897 EstimatedMass = StringMass->minMassQDiQSt
898 }
899 return EstimatedMass;
900 }
901
902 //============================================
903
904 G4double G4FTFParameters::GetProcProb( const G
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( -ProcPa
912 ProcParams[ProcN][2] * G4Exp( -ProcPa
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