| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4FTFParameters.cc,v 1.15 2010/11/15 10:02:38 vuzhinsk Exp $
>> 28 // GEANT4 tag $Name: geant4-09-04 $
27 // 29 //
28 30
29 #include <utility> <<
30 <<
31 #include "G4FTFParameters.hh" 31 #include "G4FTFParameters.hh"
32 32
33 #include "G4ios.hh" 33 #include "G4ios.hh"
34 #include "G4PhysicalConstants.hh" << 34 #include <utility>
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 35
63 G4FTFParameters::G4FTFParameters() << 36 G4FTFParameters::G4FTFParameters()
64 { << 37 {;}
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 38
94 //============================================ <<
95 39
96 void G4FTFParameters::InitForInteraction( cons << 40 G4FTFParameters::~G4FTFParameters()
97 G4in << 41 {;}
>> 42 //**********************************************************************************************
>> 43
>> 44 //G4FTFParameters::G4FTFParameters(const G4ParticleDefinition * particle,
>> 45 // G4double theA,
>> 46 // G4double theZ,
>> 47 // G4double s)
>> 48 G4FTFParameters::G4FTFParameters(const G4ParticleDefinition * particle,
>> 49 G4int theA,
>> 50 G4int theZ,
>> 51 G4double s)
98 { 52 {
99 Reset(); << 53 G4int PDGcode = particle->GetPDGEncoding();
100 << 54 G4int absPDGcode = std::abs(PDGcode);
101 G4int ProjectilePDGcode = particle->Ge << 55 G4double ProjectileMass = particle->GetPDGMass();
102 G4int ProjectileabsPDGcode = std::abs( Pr << 56 G4double TargetMass = G4Proton::Proton()->GetPDGMass();
103 G4double ProjectileMass = particle->Ge << 57
104 G4double ProjectileMass2 = ProjectileMa << 58 G4double Elab = (s - ProjectileMass*ProjectileMass - TargetMass*TargetMass)/
105 << 59 (2*TargetMass);
106 G4int ProjectileBaryonNumber( 0 ), AbsProjec << 60 G4double Plab = std::sqrt(Elab * Elab - ProjectileMass*ProjectileMass);
107 G4bool ProjectileIsNucleus = false; << 61
108 << 62 G4double Ylab=0.5*std::log((Elab+Plab)/(Elab-Plab));
109 if ( std::abs( particle->GetBaryonNumber() ) << 63
110 ProjectileIsNucleus = true; << 64 Plab/=GeV; // Uzhi 8.07.10
111 ProjectileBaryonNumber = particle->GetB << 65 G4double LogPlab = std::log( Plab );
112 AbsProjectileBaryonNumber = std::abs( Proj << 66 G4double sqrLogPlab = LogPlab * LogPlab;
113 AbsProjectileCharge = std::abs( G4in << 67
114 if ( ProjectileBaryonNumber > 1 ) { << 68 G4int NumberOfTargetProtons = theZ;
115 ProjectilePDGcode = 2212; ProjectileabsP << 69 G4int NumberOfTargetNeutrons = theA-theZ;
116 } else { << 70 // G4int NumberOfTargetProtons = (G4int) theZ;
117 ProjectilePDGcode = -2212; Projectileabs << 71 // G4int NumberOfTargetNeutrons = (G4int) theA- (G4int) theZ;
118 } << 72 G4int NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
119 ProjectileMass = G4Proton::Proton()->GetP << 73
120 ProjectileMass2 = sqr( ProjectileMass ); << 74 G4double Xtotal, Xelastic;
121 } << 75
122 << 76 if( PDGcode > 1000 ) //------Projectile is baryon --------
123 G4double TargetMass = G4Proton::Proton()->G << 77 {
124 G4double TargetMass2 = TargetMass * TargetMa << 78 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
125 << 79 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
126 G4double Plab = PlabPerParticle; << 80
127 G4double Elab = std::sqrt( Plab*Plab + Proje << 81 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
128 G4double KineticEnergy = Elab - ProjectileMa << 82 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
129 << 83
130 G4double S = ProjectileMass2 + TargetMass2 + << 84 Xtotal = ( NumberOfTargetProtons * XtotPP +
131 << 85 NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
132 #ifdef debugFTFparams << 86 Xelastic = ( NumberOfTargetProtons * XelPP +
133 G4cout << "--------- FTF Parameters -------- << 87 NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
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 } 88 }
279 << 89 else if( PDGcode < -1000 ) //------Projectile is anti_baryon --------
280 X_c = 2.0*FlowF*sqr( ProjectileMass + Ta << 90 {
281 << 91 G4double XtotPP = 38.4 + 77.6*std::pow(Plab,-0.64) + 0.26*sqrLogPlab - 1.2*LogPlab;
282 X_d = 23.3/ECMSsqr; // mb anti-quark-qu << 92 G4double XtotPN = 0. + 133.6*std::pow(Plab,-0.70) + 1.22*sqrLogPlab +13.7*LogPlab;
283 } << 93
284 << 94 G4double XelPP = 10.2 + 52.7*std::pow(Plab,-1.16) + 0.125*sqrLogPlab - 1.28*LogPlab;
285 G4double Xann_on_P( 0.0), Xann_on_N( 0.0 ) << 95 G4double XelPN = 36.5 + 0. *std::pow(Plab, 0. ) + 0. *sqrLogPlab -11.9 *LogPlab;
286 << 96
287 if ( ProjectilePDGcode == -2212 ) { << 97 Xtotal = ( NumberOfTargetProtons * XtotPP +
288 Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_ << 98 NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
289 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_ << 99 Xelastic = ( NumberOfTargetProtons * XelPP +
290 } else if ( ProjectilePDGcode == -2112 ) { << 100 NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
291 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_ << 101 }
292 Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_ << 102 else if( PDGcode == 211 ) //------Projectile is PionPlus -------
293 } else if ( ProjectilePDGcode == -3122 ) { << 103 {
294 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_ << 104 G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
295 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_ << 105 G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
296 } else if ( ProjectilePDGcode == -3112 ) { << 106
297 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_ << 107 G4double XelPiP = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
298 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_ << 108 G4double XelPiN = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
299 } else if ( ProjectilePDGcode == -3212 ) { << 109
300 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_ << 110 Xtotal = ( NumberOfTargetProtons * XtotPiP +
301 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_ << 111 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
302 } else if ( ProjectilePDGcode == -3222 ) { << 112 Xelastic = ( NumberOfTargetProtons * XelPiP +
303 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_ << 113 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
304 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_ << 114 }
305 } else if ( ProjectilePDGcode == -3312 ) { << 115 else if( PDGcode == -211 ) //------Projectile is PionMinus -------
306 Xann_on_P = X_a + X_b*1.0 + X_c*1.0 + X_ << 116 {
307 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_ << 117 G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
308 } else if ( ProjectilePDGcode == -3322 ) { << 118 G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
309 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_ << 119
310 Xann_on_N = X_a + X_b*1.0 + X_c*1.0 + X_ << 120 G4double XelPiP = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
311 } else if ( ProjectilePDGcode == -3334 ) { << 121 G4double XelPiN = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
312 Xann_on_P = X_a + X_b*0.0 + X_c*0.0 + X_ << 122
313 Xann_on_N = X_a + X_b*0.0 + X_c*0.0 + X_ << 123 Xtotal = ( NumberOfTargetProtons * XtotPiP +
314 } else { << 124 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
315 G4cout << "Unknown anti-baryon for FTF a << 125 Xelastic = ( NumberOfTargetProtons * XelPiP +
316 } << 126 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
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 } 127 }
645 128
646 } else { // The projectile is a meson as << 129 else if( PDGcode == 111 ) //------Projectile is PionZero -------
647 // Proc# A1 B1 << 130 {
648 SetParams( 0, 60.0 , 2.5 , << 131 G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
649 SetParams( 1, 6.0 , 1.0 , - << 132 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
650 SetParams( 2, 2.76, 1.2 , - << 133
651 SetParams( 3, 1.09, 0.5 , << 134 G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
652 SetParams( 4, 1.0, 0.0 , << 135 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
653 } << 136
654 << 137 G4double XelPiP =( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
655 if ( AbsProjectileBaryonNumber > 10 || N << 138 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
656 SetParams( 2, 0.0 , 0.0 , << 139 G4double XelPiN =( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
657 SetParams( 3, 0.0 , 0.0 , << 140 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
658 } << 141
659 << 142 Xtotal = ( NumberOfTargetProtons * XtotPiP +
660 SetDeltaProbAtQuarkExchange( 0.0 ); << 143 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
661 SetProbOfSameQuarkExchange( 0.0 ); << 144 Xelastic = ( NumberOfTargetProtons * XelPiP +
662 << 145 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
663 SetProjMinDiffMass( GetMinMass(particle << 146 }
664 SetProjMinNonDiffMass( GetMinMass(particle << 147 else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
665 << 148 {
666 const G4ParticleDefinition* Neutron = G4Ne << 149 G4double XtotKP = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
667 SetTarMinDiffMass( GetMinMass(Neutron)/ << 150 G4double XtotKN = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
668 SetTarMinNonDiffMass( GetMinMass(Neutron)/ << 151
669 << 152 G4double XelKP = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
670 SetAveragePt2( 0.3 ); // GeV^2 << 153 G4double XelKN = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
671 SetProbLogDistrPrD( 0.55 ); << 154
672 SetProbLogDistr( 0.55 ); << 155 Xtotal = ( NumberOfTargetProtons * XtotKP +
673 << 156 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
674 } << 157 Xelastic = ( NumberOfTargetProtons * XelKP +
675 << 158 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
676 #ifdef debugFTFparams << 159 }
677 G4cout<<"DeltaProbAtQuarkExchange "<< GetDel << 160 else if( PDGcode ==-321 ) //------Projectile is KaonMinus ------
678 G4cout<<"ProbOfSameQuarkExchange "<< GetPro << 161 {
679 G4cout<<"ProjMinDiffMass "<< GetPro << 162 G4double XtotKP = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab;
680 G4cout<<"ProjMinNonDiffMass "<< GetPro << 163 G4double XtotKN = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab;
681 G4cout<<"TarMinDiffMass "<< GetTar << 164
682 G4cout<<"TarMinNonDiffMass "<< GetTar << 165 G4double XelKP = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
683 G4cout<<"AveragePt2 "<< GetAve << 166 G4double XelKN = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
684 G4cout<<"ProbLogDistrPrD "<< GetPro << 167
685 G4cout<<"ProbLogDistrTrD "<< GetPro << 168 Xtotal = ( NumberOfTargetProtons * XtotKP +
686 #endif << 169 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
687 << 170 Xelastic = ( NumberOfTargetProtons * XelKP +
688 // Set parameters of nuclear destruction << 171 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
689 << 172 }
690 if ( ProjectileabsPDGcode < 1000 ) { // Mes << 173 else if((PDGcode == 311) || (PDGcode == 130) || (PDGcode == 310))//Projectile is KaonZero
691 << 174 {
692 const G4int indexTune = G4FTFTunings::Inst << 175 G4double XtotKP =( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
693 << 176 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
694 SetMaxNumberOfCollisions( Plab, 2.0 ); // << 177 G4double XtotKN =( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
695 // << 178 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
696 // target destruction << 179
697 // << 180 G4double XelKP =( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
698 /* original code ---> << 181 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
699 SetCofNuclearDestruction( 0.00481*G4double << 182 G4double XelKN =( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
700 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + << 183 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
701 << 184 Xtotal = ( NumberOfTargetProtons * XtotKP +
702 SetR2ofNuclearDestruction( 1.5*fermi*fermi << 185 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
703 SetDofNuclearDestruction( 0.3 ); << 186 Xelastic = ( NumberOfTargetProtons * XelKP +
704 SetPt2ofNuclearDestruction( ( 0.035 + 0.04 << 187 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
705 ( 1.0 << 188 }
706 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV << 189 else //------Projectile is undefined, Nucleon assumed
707 SetExcitationEnergyPerWoundedNucleon( 40.0 << 190 {
708 */ << 191 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
709 double coeff = fArrayParCollMesonProj[inde << 192 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
710 // << 193
711 // NOTE (JVY): Set this switch to false/tr << 194 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
712 // << 195 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
713 if ( fArrayParCollMesonProj[indexTune].IsN << 196
714 { << 197 Xtotal = ( NumberOfTargetProtons * XtotPP +
715 coeff *= G4double(NumberOfTargetNucleons << 198 NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
716 } << 199 Xelastic = ( NumberOfTargetProtons * XelPP +
717 double exfactor = G4Exp( fArrayParCollMeso << 200 NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
718 * (Ylab-fArrayParCo << 201 };
719 coeff *= exfactor; << 202
720 coeff /= ( 1.+ exfactor ); << 203 // Xtotal and Xelastic in mb
721 << 204
722 SetCofNuclearDestruction( coeff ); << 205 // For Pi- P interactions only!
723 << 206 if(std::abs(Plab-1.4) < 0.05) {Xtotal=3.500599e+01; Xelastic= 1.150032e+01;}
724 SetR2ofNuclearDestruction( fArrayParCollMe << 207 if(std::abs(Plab-1.5) < 0.05) {Xtotal=3.450591e+01; Xelastic= 1.050038e+01;}
725 SetDofNuclearDestruction( fArrayParCollMes << 208 if(std::abs(Plab-1.6) < 0.05) {Xtotal=3.430576e+01; Xelastic= 9.800433e+00;}
726 coeff = fArrayParCollMesonProj[indexTune]. << 209 if(std::abs(Plab-1.7) < 0.05) {Xtotal=3.455560e+01; Xelastic= 9.300436e+00;}
727 exfactor = G4Exp( fArrayParCollMesonProj[i << 210 if(std::abs(Plab-1.8) < 0.05) {Xtotal=3.480545e+01; Xelastic= 8.800438e+00;}
728 * (Ylab-fArrayParCollMeson << 211 if(std::abs(Plab-2.0) < 0.05) {Xtotal=3.570503e+01; Xelastic= 8.200370e+00;}
729 coeff *= exfactor; << 212 if(std::abs(Plab-2.2) < 0.05) {Xtotal=3.530495e+01; Xelastic= 7.800362e+00;}
730 coeff /= ( 1. + exfactor ); << 213 if(std::abs(Plab-2.5) < 0.05) {Xtotal=3.410484e+01; Xelastic= 7.350320e+00;}
731 SetPt2ofNuclearDestruction( (fArrayParColl << 214 if(std::abs(Plab-2.75) < 0.05){Xtotal=3.280479e+01; Xelastic= 7.050273e+00;}
732 << 215 if(std::abs(Plab-3.0) < 0.05) {Xtotal=3.180473e+01; Xelastic= 6.800258e+00;}
733 SetMaxPt2ofNuclearDestruction( fArrayParCo << 216 if(std::abs(Plab-4.0) < 0.05) {Xtotal=2.910441e+01; Xelastic= 6.100229e+00;}
734 SetExcitationEnergyPerWoundedNucleon( fArr << 217 if(std::abs(Plab-5.0) < 0.05) {Xtotal=2.820372e+01; Xelastic= 5.700275e+00;}
735 << 218 if(std::abs(Plab-6.0) < 0.05) {Xtotal=2.760367e+01; Xelastic= 5.400255e+00;}
736 } else if ( ProjectilePDGcode == -2212 || << 219 if(std::abs(Plab-7.0) < 0.05) {Xtotal=2.725366e+01; Xelastic= 5.150256e+00;}
737 << 220 if(std::abs(Plab-8.0) < 0.05) {Xtotal=2.690365e+01; Xelastic= 4.900258e+00;}
738 SetMaxNumberOfCollisions( Plab, 2.0 ); << 221 if(std::abs(Plab-10.0) < 0.05){Xtotal=2.660342e+01; Xelastic= 4.600237e+00;}
739 << 222 if(std::abs(Plab-12.0) < 0.05){Xtotal=2.632341e+01; Xelastic= 4.480229e+00;}
740 SetCofNuclearDestruction( 0.00481*G4double << 223 if(std::abs(Plab-14.0) < 0.05){Xtotal=2.604340e+01; Xelastic= 4.360221e+00;}
741 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G << 224 if(std::abs(Plab-20.0) < 0.05){Xtotal=2.520337e+01; Xelastic= 4.000197e+00;}
742 SetR2ofNuclearDestruction( 1.5*fermi*fermi << 225 if(std::abs(Plab-30.0) < 0.05){Xtotal=2.505334e+01; Xelastic= 3.912679e+00;}
743 SetDofNuclearDestruction( 0.3 ); << 226 //
744 SetPt2ofNuclearDestruction( ( 0.035 + 0.04 << 227 //----------- Geometrical parameters ------------------------------------------------
745 ( 1.0 << 228 SetTotalCrossSection(Xtotal);
746 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV << 229 SetElastisCrossSection(Xelastic);
747 SetExcitationEnergyPerWoundedNucleon( 40.0 << 230 SetInelasticCrossSection(Xtotal-Xelastic);
748 if ( Plab < 2.0 ) { // 2 GeV/c << 231
749 // For slow anti-baryon we have to garan << 232 //G4cout<<"Plab Xtotal, Xelastic Xinel "<<Plab<<" "<<Xtotal<<" "<<Xelastic<<Xtotal-Xelastic)<<G4endl;
750 SetCofNuclearDestruction( 0.0 ); << 233 // // Interactions with elastic and inelastic collisions
751 SetR2ofNuclearDestruction( 1.5*fermi*fer << 234 SetProbabilityOfElasticScatt(Xtotal, Xelastic);
752 << 235 SetRadiusOfHNinteractions2(Xtotal/pi/10.);
753 SetDofNuclearDestruction( 0.01 ); << 236 //
754 SetPt2ofNuclearDestruction( 0.035*GeV*Ge << 237 /* //==== No elastic scattering ============================
755 SetMaxPt2ofNuclearDestruction( 0.04*GeV* << 238 SetProbabilityOfElasticScatt(Xtotal, 0.);
756 } << 239 SetRadiusOfHNinteractions2((Xtotal-Xelastic)/pi/10.);
757 << 240 */ //=======================================================
758 } else { // Projectile baryon assumed << 241
759 << 242 //-----------------------------------------------------------------------------------
760 // Below, in the call to the G4FTFTunings: << 243
761 // as projectile, instead of the real one, << 244 SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
762 // destruction, we assume for this categor << 245 // (GeV/c)^(-2))
763 // as for "baryon". << 246 //-----------------------------------------------------------------------------------
764 const G4int indexTune = G4FTFTunings::Inst << 247 SetGamma0( GetSlope()*Xtotal/10./2./pi );
765 << 248
766 // NOTE (JVY) FIXME !!! Will decide later << 249 //----------- Parameters of elastic scattering --------------------------------------
767 // << 250 // Gaussian parametrization of
768 SetMaxNumberOfCollisions( Plab, 2.0 ); << 251 // elastic scattering amplitude assumed
769 << 252 SetAvaragePt2ofElasticScattering(1./(Xtotal*Xtotal/16./pi/Xelastic/0.3894)*GeV*GeV);
770 // projectile destruction - does NOT reall << 253
771 // << 254 //----------- Parameters of excitations ---------------------------------------------
772 double coeff = 0.; << 255 if( PDGcode > 1000 ) //------Projectile is baryon --------
773 coeff = fArrayParCollBaryonProj[indexTune] << 256 {
774 // << 257 SetMagQuarkExchange(1.84);//(3.63);
775 // NOTE (JVY): Set this switch to false/tr << 258 SetSlopeQuarkExchange(0.7);//(1.2);
776 // << 259 SetDeltaProbAtQuarkExchange(0.);
777 if ( fArrayParCollBaryonProj[indexTune].Is << 260 if(NumberOfTargetNucleons > 26) {SetProbOfSameQuarkExchange(1.);}
778 { << 261 else {SetProbOfSameQuarkExchange(0.);}
779 coeff *= G4double(AbsProjectileBaryonNum << 262
780 } << 263 SetProjMinDiffMass(1.16); // GeV
781 double exfactor = G4Exp( fArrayParCollBary << 264 SetProjMinNonDiffMass(1.16); // GeV
782 (Ylab-fArrayParCollBaryonProj[index << 265 SetProbabilityOfProjDiff(0.805*std::exp(-0.35*Ylab));// 0.5
783 coeff *= exfactor; << 266
784 coeff /= ( 1.+ exfactor ); << 267 SetTarMinDiffMass(1.16); // GeV
785 SetCofNuclearDestructionPr( coeff ); << 268 SetTarMinNonDiffMass(1.16); // GeV
786 << 269 SetProbabilityOfTarDiff(0.805*std::exp(-0.35*Ylab));// 0.5
787 // target desctruction << 270
788 // << 271 SetAveragePt2(0.15); // 0.15 GeV^2
789 coeff = fArrayParCollBaryonProj[indexTune] << 272 }
790 // << 273 if( PDGcode < -1000 ) //------Projectile is anti_baryon --------
791 // NOTE (JVY): Set this switch to false/tr << 274 {
792 // << 275 SetMagQuarkExchange(0.);
793 if ( fArrayParCollBaryonProj[indexTune].Is << 276 SetSlopeQuarkExchange(0.);
794 { << 277 SetDeltaProbAtQuarkExchange(0.);
795 coeff *= G4double(NumberOfTargetNucleons << 278 SetProbOfSameQuarkExchange(0.);
796 } << 279
797 exfactor = G4Exp( fArrayParCollBaryonProj[ << 280 SetProjMinDiffMass(1.16); // GeV
798 (Ylab-fArrayParCollBaryonProj[indexT << 281 SetProjMinNonDiffMass(1.16); // GeV
799 coeff *= exfactor; << 282 SetProbabilityOfProjDiff(0.805*std::exp(-0.35*Ylab));// 0.5
800 coeff /= ( 1.+ exfactor ); << 283
801 SetCofNuclearDestruction( coeff ); << 284 SetTarMinDiffMass(1.16); // GeV
802 << 285 SetTarMinNonDiffMass(1.16); // GeV
803 SetR2ofNuclearDestruction( fArrayParCollBa << 286 SetProbabilityOfTarDiff(0.805*std::exp(-0.35*Ylab));// 0.5
804 SetDofNuclearDestruction( fArrayParCollBar << 287
805 << 288 SetAveragePt2(0.15); // 0.15 GeV^2
806 coeff = fArrayParCollBaryonProj[indexTune] << 289 }
807 exfactor = G4Exp( fArrayParCollBaryonProj[ << 290 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion -----------
808 (Ylab-fArrayParCollBaryonProj[indexT << 291 {
809 coeff *= exfactor; << 292 SetMagQuarkExchange(240.);
810 coeff /= ( 1. + exfactor ); << 293 SetSlopeQuarkExchange(2.);
811 SetPt2ofNuclearDestruction( (fArrayParColl << 294 SetDeltaProbAtQuarkExchange(0.56); //(0.35);
812 << 295
813 SetMaxPt2ofNuclearDestruction( fArrayParCo << 296 SetProjMinDiffMass(0.5); // GeV
814 SetExcitationEnergyPerWoundedNucleon( fArr << 297 SetProjMinNonDiffMass(0.5); // GeV 0.3
815 << 298 SetProbabilityOfProjDiff(0.);//(0.*0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
816 } << 299
817 << 300 SetTarMinDiffMass(1.16); // GeV
818 #ifdef debugFTFparams << 301 SetTarMinNonDiffMass(1.16); // GeV
819 G4cout<<"CofNuclearDestructionPr " << 302 // SetProbabilityOfTarDiff(1.);//(2.*0.62*std::pow(s/GeV/GeV,-0.51));
820 G4cout<<"CofNuclearDestructionTr " << 303 // SetProbabilityOfTarDiff(2.6*std::exp(-0.46*Ylab));
821 G4cout<<"R2ofNuclearDestruction " << 304 SetProbabilityOfTarDiff(0.8*std::exp(-0.6*(Ylab-3.)));
822 G4cout<<"DofNuclearDestruction " << 305
823 G4cout<<"Pt2ofNuclearDestruction " << 306 SetAveragePt2(0.3); // GeV^2
824 G4cout<<"ExcitationEnergyPerWoundedNucleon " << 307 }
825 #endif << 308 else if( (absPDGcode == 321) || (PDGcode == 311) ||
826 << 309 (PDGcode == 130) || (PDGcode == 310)) //Projectile is Kaon
827 //SetCofNuclearDestruction( 0.47*G4Exp( 2.0* << 310 {
828 //SetPt2ofNuclearDestruction( ( 0.035 + 0.1* << 311 // Must be corrected, taken from PiN
829 << 312 SetMagQuarkExchange(120.);
830 //SetMagQuarkExchange( 120.0 ); // 210 << 313 SetSlopeQuarkExchange(2.0);
831 //SetSlopeQuarkExchange( 2.0 ); << 314 SetDeltaProbAtQuarkExchange(0.6);
832 //SetDeltaProbAtQuarkExchange( 0.6 ); << 315
833 //SetProjMinDiffMass( 0.7 ); // GeV << 316 SetProjMinDiffMass(0.7); // GeV 1.1
834 //SetProjMinNonDiffMass( 0.7 ); // GeV << 317 SetProjMinNonDiffMass(0.7); // GeV
835 //SetProbabilityOfProjDiff( 0.0); // 0.8 << 318 SetProbabilityOfProjDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
836 //SetTarMinDiffMass( 1.1 ); // GeV << 319
837 //SetTarMinNonDiffMass( 1.1 ); // GeV << 320 SetTarMinDiffMass(1.1); // GeV
838 //SetProbabilityOfTarDiff( 0.0 ); // 0.8 << 321 SetTarMinNonDiffMass(1.1); // GeV
839 << 322 SetProbabilityOfTarDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
840 //SetAveragePt2( 0.0 ); // GeV << 323
841 //------------------------------------ << 324 SetAveragePt2(0.3); // GeV^2
842 //SetProbabilityOfElasticScatt( 1.0, 1.0); << 325 }
843 //SetProbabilityOfProjDiff( 1.0*0.62*G4Pow:: << 326 else //------Projectile is undefined,
844 //SetProbabilityOfTarDiff( 4.0*0.62*G4Pow::G << 327 //------Nucleon assumed
845 //SetAveragePt2( 0.3 ); << 328 {
846 //SetAvaragePt2ofElasticScattering( 0.0 ); << 329 SetMagQuarkExchange(3.5);
847 << 330 SetSlopeQuarkExchange(1.0);
848 //SetMaxNumberOfCollisions( Plab, 6.0 ); //( << 331 SetDeltaProbAtQuarkExchange(0.1);
849 //SetAveragePt2( 0.15 ); << 332
850 //SetCofNuclearDestruction(-1.);//( 0.75 ); << 333 SetProjMinDiffMass((particle->GetPDGMass()+160.*MeV)/GeV);
851 //SetExcitationEnergyPerWoundedNucleon(0.);/ << 334 SetProjMinNonDiffMass((particle->GetPDGMass()+160.*MeV)/GeV);
852 //SetDofNuclearDestruction(0.);//( 0.2 ); // << 335 SetProbabilityOfProjDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
853 << 336
854 /* << 337 SetTarMinDiffMass(1.1); // GeV
855 SetAveragePt2(0.3); << 338 SetTarMinNonDiffMass(1.1); // GeV
856 SetCofNuclearDestructionPr(0.); << 339 SetProbabilityOfTarDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
857 SetCofNuclearDestruction(0.); // << 340
858 SetExcitationEnergyPerWoundedNucleon(0.); // << 341 SetAveragePt2(0.3); // GeV^2
859 SetDofNuclearDestruction(0.); // << 342 }
860 SetPt2ofNuclearDestruction(0.); // << 343
861 */ << 344 // ---------- Set parameters of a string kink -------------------------------
862 << 345 SetPt2Kink(6.*GeV*GeV);
863 //SetExcitationEnergyPerWoundedNucleon(0.001 << 346 G4double Puubar(1./3.), Pddbar(1./3.), Pssbar(1./3.); // SU(3) symmetry
864 //SetPt2Kink( 0.0*GeV*GeV ); << 347 // G4double Puubar(0.41 ), Pddbar(0.41 ), Pssbar(0.18 ); // Broken SU(3) symmetry
865 << 348 SetQuarkProbabilitiesAtGluonSplitUp(Puubar, Pddbar, Pssbar);
866 //SetRadiusOfHNinteractions2( Xtotal/pi/10.0 << 349
867 //SetRadiusOfHNinteractions2( (Xtotal - Xela << 350 // --------- Set parameters of nuclear destruction--------------------
868 //SetProbabilityOfElasticScatt( 1.0, 0.0); << 351
869 /* << 352 if( absPDGcode < 1000 )
870 G4cout << "Pt2 " << GetAveragePt2()<<" "<<Ge << 353 {
871 G4cout << "Cnd " << GetCofNuclearDestruction << 354 SetMaxNumberOfCollisions(1000.,1.); //(Plab,2.); //3.); ##############################
872 G4cout << "Dnd " << GetDofNuclearDestruction << 355
873 G4cout << "Pt2 " << GetPt2ofNuclearDestructi << 356 // SetCofNuclearDestruction(0.); //1.0); // for meson projectile
874 */ << 357 // SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/(1.+std::exp(4.*(Ylab-2.1))));
875 << 358 //G4cout<<Ylab<<" "<<0.62*std::exp(4.*(Ylab-4.5))/(1.+std::exp(4.*(Ylab-4.5)))<<G4endl;
876 } << 359 //G4int Uzhi; G4cin>>Uzhi;
877 << 360
878 //============================================ << 361 // SetMaxNumberOfCollisions(Plab,2.); //4.); // ##############################
879 << 362
880 G4double G4FTFParameters::GetMinMass( const G4 << 363 SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/(1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
881 // The code is used for estimating the mini << 364 //------------------------------------------
882 // The indices used for minMassQDiQStr must << 365 // SetDofNuclearDestruction(0.4);
883 // quarks: d, u, s, c and b; enforcing this << 366 // SetPt2ofNuclearDestruction(0.17*GeV*GeV);
884 G4double EstimatedMass = 0.0; << 367 // SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
885 G4int partID = std::abs(aParticle->GetPDGEn << 368
886 G4int Qleft = std::max( partID/100, 1 << 369 // SetExcitationEnergyPerWoundedNucleon(100*MeV);
887 G4int Qright = std::max( (partID/ 10)%10, 1 << 370 SetDofNuclearDestruction(0.4);
888 if ( Qleft < 6 && Qright < 6 ) { << 371 SetPt2ofNuclearDestruction((0.035+
889 EstimatedMass = StringMass->minMassQQbarS << 372 0.04*std::exp(4.*(Ylab-2.5))/(1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
890 } else if ( Qleft < 6 && Qright > 6 ) { << 373 SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
891 G4int q1 = std::max( std::min( Qright/10, << 374
892 G4int q2 = std::max( std::min( Qright%10, << 375 SetExcitationEnergyPerWoundedNucleon(75.*MeV);
893 EstimatedMass = StringMass->minMassQDiQSt << 376 } else // for baryon projectile
894 } else if ( Qleft > 6 && Qright < 6 ) { << 377 {
895 G4int q1 = std::max( std::min( Qleft/10, << 378 SetMaxNumberOfCollisions(Plab,2.); //4.); // ##############################
896 G4int q2 = std::max( std::min( Qleft%10, << 379
897 EstimatedMass = StringMass->minMassQDiQSt << 380 SetCofNuclearDestruction(1.*std::exp(4.*(Ylab-2.1))/(1.+std::exp(4.*(Ylab-2.1)))); //0.62 1.0
898 } << 381 //G4cout<<Ylab<<" "<<0.62*std::exp(4.*(Ylab-2.1))/(1.+std::exp(4.*(Ylab-2.1)))<<G4endl;
899 return EstimatedMass; << 382 //G4int Uzhi; G4cin>>Uzhi;
900 } << 383
901 << 384 SetDofNuclearDestruction(0.4);
902 //============================================ << 385 SetPt2ofNuclearDestruction((0.035+
903 << 386 0.04*std::exp(4.*(Ylab-2.5))/(1.+std::exp(4.*(Ylab-2.5))))*GeV*GeV); //0.09
904 G4double G4FTFParameters::GetProcProb( const G << 387 SetMaxPt2ofNuclearDestruction(1.0*GeV*GeV);
905 G4double Prob( 0.0 ); << 388
906 if ( y < ProcParams[ProcN][6] ) { << 389 SetExcitationEnergyPerWoundedNucleon(75.*MeV);
907 Prob = ProcParams[ProcN][5]; << 390 }
908 if (Prob < 0.) Prob=0.; << 391
909 return Prob; << 392 SetR2ofNuclearDestruction(1.5*fermi*fermi);
910 } << 393
911 Prob = ProcParams[ProcN][0] * G4Exp( -ProcPa << 394 //SetCofNuclearDestruction(0.47*std::exp(2.*(Ylab-2.5))/(1.+std::exp(2.*(Ylab-2.5))));
912 ProcParams[ProcN][2] * G4Exp( -ProcPa << 395 //SetPt2ofNuclearDestruction((0.035+0.1*std::exp(4.*(Ylab-3.))/(1.+std::exp(4.*(Ylab-3.))))*GeV*GeV);
913 ProcParams[ProcN][4]; << 396
914 if (Prob < 0.) Prob=0.; << 397 //SetMagQuarkExchange(120.); // 210. PipP
915 return Prob; << 398 //SetSlopeQuarkExchange(2.0);
916 } << 399 //SetDeltaProbAtQuarkExchange(0.6);
917 << 400 //SetProjMinDiffMass(0.7); // GeV 1.1
918 //============================================ << 401 //SetProjMinNonDiffMass(0.7); // GeV
919 << 402 //SetProbabilityOfProjDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
920 G4FTFParameters::~G4FTFParameters() { << 403 //SetTarMinDiffMass(1.1); // GeV
921 if ( StringMass ) delete StringMass; << 404 //SetTarMinNonDiffMass(1.1); // GeV
922 } << 405 //SetProbabilityOfTarDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
923 << 406 //
924 //============================================ << 407 //SetAveragePt2(0.3); // GeV^2
925 << 408 //------------------------------------
926 void G4FTFParameters::Reset() << 409 //SetProbabilityOfElasticScatt(1.,1.); //(Xtotal, Xelastic);
927 { << 410 //SetProbabilityOfProjDiff(1.*0.62*std::pow(s/GeV/GeV,-0.51)); // 0->1
928 FTFhNcmsEnergy = 0.0; << 411 //SetProbabilityOfTarDiff(4.*0.62*std::pow(s/GeV/GeV,-0.51)); // 2->4
929 FTFXtotal = 0.0; << 412 //SetAveragePt2(0.3); //(0.15);
930 FTFXelastic = 0.0; << 413 //SetAvaragePt2ofElasticScattering(0.);
931 FTFXinelastic = 0.0; << 414
932 FTFXannihilation = 0.0; << 415 //SetMaxNumberOfCollisions(4.*(Plab+0.01),Plab); //6.); // ##############################
933 ProbabilityOfAnnihilation = 0.0; << 416 //SetCofNuclearDestruction(0.2); //(0.4);
934 ProbabilityOfElasticScatt = 0.0; << 417 //SetExcitationEnergyPerWoundedNucleon(0.*MeV); //(75.*MeV);
935 RadiusOfHNinteractions2 = 0.0; << 418 //SetDofNuclearDestruction(0.4); //(0.4);
936 FTFSlope = 0.0; << 419 //SetPt2ofNuclearDestruction(0.1*GeV*GeV); //(0.168*GeV*GeV);
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 420
>> 421 }
>> 422 //**********************************************************************************************
968 423