| 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.4 2008/12/18 13:02:00 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-09-02-patch-03 $
27 // 29 //
28 30
29 #include <utility> <<
30 <<
31 #include "G4FTFParameters.hh" 31 #include "G4FTFParameters.hh"
32 32
33 #include "G4ios.hh" << 33 G4FTFParameters::G4FTFParameters()
34 #include "G4PhysicalConstants.hh" << 34 {;}
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 35
282 X_d = 23.3/ECMSsqr; // mb anti-quark-qu <<
283 } <<
284 36
285 G4double Xann_on_P( 0.0), Xann_on_N( 0.0 ) << 37 G4FTFParameters::~G4FTFParameters()
>> 38 {;}
>> 39 //**********************************************************************************************
>> 40
>> 41 G4FTFParameters::G4FTFParameters(const G4ParticleDefinition * particle,
>> 42 G4double theA,
>> 43 G4double theZ,
>> 44 G4double s)
>> 45 {
>> 46 G4int PDGcode = particle->GetPDGEncoding();
>> 47 G4int absPDGcode = std::abs(PDGcode);
>> 48 G4double Elab = (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
>> 49 G4double Plab = std::sqrt(Elab * Elab - 0.88);
>> 50
>> 51 G4double LogPlab = std::log( Plab );
>> 52 G4double sqrLogPlab = LogPlab * LogPlab;
>> 53
>> 54 //G4cout<<"G4FTFParameters Plab "<<Plab<<G4endl;
>> 55
>> 56 G4int NumberOfTargetProtons = (G4int) theZ;
>> 57 G4int NumberOfTargetNeutrons = (G4int) theA- (G4int) theZ;
>> 58 G4int NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
>> 59
>> 60 G4double Xtotal, Xelastic;
>> 61
>> 62 if( absPDGcode > 1000 ) //------Projectile is baryon --------
>> 63 {
>> 64 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
>> 65 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
>> 66
>> 67 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
>> 68 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
>> 69
>> 70 Xtotal = ( NumberOfTargetProtons * XtotPP +
>> 71 NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
>> 72 Xelastic = ( NumberOfTargetProtons * XelPP +
>> 73 NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
>> 74 }
>> 75 else if( PDGcode == 211 ) //------Projectile is PionPlus -------
>> 76 {
>> 77 G4double XtotPiP = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
>> 78 G4double XtotPiN = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
>> 79
>> 80 G4double XelPiP = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
>> 81 G4double XelPiN = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
>> 82
>> 83 Xtotal = ( NumberOfTargetProtons * XtotPiP +
>> 84 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
>> 85 Xelastic = ( NumberOfTargetProtons * XelPiP +
>> 86 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
>> 87 }
>> 88 else if( PDGcode == -211 ) //------Projectile is PionMinus -------
>> 89 {
>> 90 G4double XtotPiP = 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab;
>> 91 G4double XtotPiN = 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab;
>> 92
>> 93 G4double XelPiP = 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab;
>> 94 G4double XelPiN = 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab;
>> 95
>> 96 Xtotal = ( NumberOfTargetProtons * XtotPiP +
>> 97 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
>> 98 Xelastic = ( NumberOfTargetProtons * XelPiP +
>> 99 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
>> 100 }
286 101
287 if ( ProjectilePDGcode == -2212 ) { << 102 else if( PDGcode == 111 ) //------Projectile is PionZero -------
288 Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_ << 103 {
289 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_ << 104 G4double XtotPiP =(16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab + //Pi+
290 } else if ( ProjectilePDGcode == -2112 ) { << 105 33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab)/2; //Pi-
291 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_ << 106
292 Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_ << 107 G4double XtotPiN =(33.0 + 14.0 *std::pow(Plab,-1.36) + 0.456*sqrLogPlab - 4.03*LogPlab + //Pi+
293 } else if ( ProjectilePDGcode == -3122 ) { << 108 16.4 + 19.3 *std::pow(Plab,-0.42) + 0.19 *sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
294 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_ << 109
295 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_ << 110 G4double XelPiP =( 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab + //Pi+
296 } else if ( ProjectilePDGcode == -3112 ) { << 111 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
297 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_ << 112 G4double XelPiN =( 1.76 + 11.2*std::pow(Plab,-0.64) + 0.043*sqrLogPlab - 0.0 *LogPlab + //Pi+
298 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_ << 113 0.0 + 11.4*std::pow(Plab,-0.40) + 0.079*sqrLogPlab - 0.0 *LogPlab)/2; //Pi-
299 } else if ( ProjectilePDGcode == -3212 ) { << 114
300 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_ << 115 Xtotal = ( NumberOfTargetProtons * XtotPiP +
301 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_ << 116 NumberOfTargetNeutrons * XtotPiN ) / NumberOfTargetNucleons;
302 } else if ( ProjectilePDGcode == -3222 ) { << 117 Xelastic = ( NumberOfTargetProtons * XelPiP +
303 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_ << 118 NumberOfTargetNeutrons * XelPiN ) / NumberOfTargetNucleons;
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 } 119 }
>> 120 else if( PDGcode == 321 ) //------Projectile is KaonPlus -------
>> 121 {
>> 122 G4double XtotKP = 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab;
>> 123 G4double XtotKN = 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab;
>> 124
>> 125 G4double XelKP = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
>> 126 G4double XelKN = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
>> 127
>> 128 Xtotal = ( NumberOfTargetProtons * XtotKP +
>> 129 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
>> 130 Xelastic = ( NumberOfTargetProtons * XelKP +
>> 131 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
>> 132 }
>> 133 else if( PDGcode ==-321 ) //------Projectile is KaonMinus ------
>> 134 {
>> 135 G4double XtotKP = 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab;
>> 136 G4double XtotKN = 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab;
>> 137
>> 138 G4double XelKP = 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab;
>> 139 G4double XelKN = 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab;
>> 140
>> 141 Xtotal = ( NumberOfTargetProtons * XtotKP +
>> 142 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
>> 143 Xelastic = ( NumberOfTargetProtons * XelKP +
>> 144 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
>> 145 }
>> 146 else if( PDGcode == 311 ) //------Projectile is KaonZero ------
>> 147 {
>> 148 G4double XtotKP =( 18.1 + 0. *std::pow(Plab, 0. ) + 0.26 *sqrLogPlab - 1.0 *LogPlab + //K+
>> 149 32.1 + 0. *std::pow(Plab, 0. ) + 0.66 *sqrLogPlab - 5.6 *LogPlab)/2; //K-
>> 150 G4double XtotKN =( 18.7 + 0. *std::pow(Plab, 0. ) + 0.21 *sqrLogPlab - 0.89*LogPlab + //K+
>> 151 25.2 + 0. *std::pow(Plab, 0. ) + 0.38 *sqrLogPlab - 2.9 *LogPlab)/2; //K-
>> 152
>> 153 G4double XelKP =( 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab + //K+
>> 154 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab)/2; //K-
>> 155 G4double XelKN =( 7.3 + 0. *std::pow(Plab,-0. ) + 0.29 *sqrLogPlab - 2.4 *LogPlab + //K+
>> 156 5.0 + 8.1*std::pow(Plab,-1.8 ) + 0.16 *sqrLogPlab - 1.3 *LogPlab)/2; //K-
>> 157 Xtotal = ( NumberOfTargetProtons * XtotKP +
>> 158 NumberOfTargetNeutrons * XtotKN ) / NumberOfTargetNucleons;
>> 159 Xelastic = ( NumberOfTargetProtons * XelKP +
>> 160 NumberOfTargetNeutrons * XelKN ) / NumberOfTargetNucleons;
>> 161 }
>> 162 else //------Projectile is undefined, Nucleon assumed
>> 163 {
>> 164 G4double XtotPP = 48.0 + 0. *std::pow(Plab, 0. ) + 0.522*sqrLogPlab - 4.51*LogPlab;
>> 165 G4double XtotPN = 47.3 + 0. *std::pow(Plab, 0. ) + 0.513*sqrLogPlab - 4.27*LogPlab;
>> 166
>> 167 G4double XelPP = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
>> 168 G4double XelPN = 11.9 + 26.9*std::pow(Plab,-1.21) + 0.169*sqrLogPlab - 1.85*LogPlab;
>> 169
>> 170 Xtotal = ( NumberOfTargetProtons * XtotPP +
>> 171 NumberOfTargetNeutrons * XtotPN ) / NumberOfTargetNucleons;
>> 172 Xelastic = ( NumberOfTargetProtons * XelPP +
>> 173 NumberOfTargetNeutrons * XelPN ) / NumberOfTargetNucleons;
>> 174 };
>> 175
>> 176 // Xtotal and Xelastic in mb
>> 177
>> 178 //----------- Geometrical parameters ------------------------------------------------
>> 179 SetTotalCrossSection(Xtotal);
>> 180 SetElastisCrossSection(Xelastic);
>> 181 SetInelasticCrossSection(Xtotal-Xelastic);
>> 182
>> 183 // // Interactions with elastic ans inelastic collisions
>> 184 SetProbabilityOfElasticScatt(Xtotal, Xelastic);
>> 185 SetRadiusOfHNinteractions2(Xtotal/pi/10.);
>> 186 //
>> 187 /* //==== No elastic scattering ============================
>> 188 SetProbabilityOfElasticScatt(Xtotal, 0.);
>> 189 SetRadiusOfHNinteractions2((Xtotal-Xelastic)/pi/10.);
>> 190 */ //=======================================================
>> 191
>> 192 //G4cout<<" Rnn "<<Xtotal/pi/10.<<" "<<Xtotal/pi/10.*fermi*fermi<<G4endl;
>> 193 //G4cout<<"G4FTFParameters Xt Xel MeV "<<Xtotal<<" "<<Xelastic<<" "<<GeV<<G4endl;
>> 194
>> 195 //-----------------------------------------------------------------------------------
>> 196 SetSlope( Xtotal*Xtotal/16./pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
>> 197 // (GeV/c)^(-2))
>> 198 //G4cout<<"G4FTFParameters Slope "<<GetSlope()<<G4endl;
>> 199 //-----------------------------------------------------------------------------------
>> 200 SetGamma0( GetSlope()*Xtotal/10./2./pi );
>> 201
>> 202 //----------- Parameters of elastic scattering --------------------------------------
>> 203 // Gaussian parametrization of
>> 204 // elastic scattering amplitude assumed
>> 205 SetAvaragePt2ofElasticScattering(1./(Xtotal*Xtotal/16./pi/Xelastic/0.3894)*GeV*GeV);
>> 206
>> 207 //----------- Parameters of excitations ---------------------------------------------
>> 208 if( absPDGcode > 1000 ) //------Projectile is baryon --------
>> 209 {
>> 210 SetProjMinDiffMass(1.1); // GeV
>> 211 SetProjMinNonDiffMass(1.1); // GeV
>> 212 SetProbabilityOfProjDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
>> 213
>> 214 SetTarMinDiffMass(1.1); // GeV
>> 215 SetTarMinNonDiffMass(1.1); // GeV
>> 216 SetProbabilityOfTarDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
>> 217
>> 218 SetAveragePt2(0.3); // GeV^2
>> 219 }
>> 220 else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion -----------
>> 221 {
>> 222 SetProjMinDiffMass(0.5); // GeV
>> 223 SetProjMinNonDiffMass(0.3); // GeV
>> 224 SetProbabilityOfProjDiff(0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
>> 225
>> 226 SetTarMinDiffMass(1.1); // GeV
>> 227 SetTarMinNonDiffMass(1.1); // GeV
>> 228 SetProbabilityOfTarDiff(0.62*std::pow(s/GeV/GeV,-0.51)); // 40/32 X-dif/X-inel
>> 229
>> 230 /*
>> 231 SetProjMinDiffMass(0.5);
>> 232 SetProjMinNonDiffMass(0.3); // Uzhi 12.06.08
>> 233 SetProbabilityOfProjDiff(0.05);
>> 234 SetProbabilityOfTarDiff(0.05);
>> 235 */
>> 236 SetAveragePt2(0.3); // GeV^2
>> 237 }
>> 238 else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
>> 239 {
>> 240 SetProjMinDiffMass(0.7); // GeV 1.1
>> 241 SetProjMinNonDiffMass(0.7); // GeV
>> 242 SetProbabilityOfProjDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
>> 243
>> 244 SetTarMinDiffMass(1.1); // GeV
>> 245 SetTarMinNonDiffMass(1.1); // GeV
>> 246 SetProbabilityOfTarDiff(0.85*std::pow(s/GeV/GeV,-0.5)); // 40/32 X-dif/X-inel
>> 247
>> 248 SetAveragePt2(0.3); // GeV^2
>> 249 }
>> 250 else //------Projectile is undefined,
>> 251 //------Nucleon assumed
>> 252 {
>> 253 SetProjMinDiffMass((particle->GetPDGMass()+160.*MeV)/GeV);
>> 254 SetProjMinNonDiffMass((particle->GetPDGMass()+160.*MeV)/GeV);
>> 255 SetProbabilityOfProjDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
>> 256
>> 257 SetTarMinDiffMass(1.1); // GeV
>> 258 SetTarMinNonDiffMass(1.1); // GeV
>> 259 SetProbabilityOfTarDiff(0.95*std::pow(s/GeV/GeV,-0.35)); // 40/32 X-dif/X-inel
>> 260
>> 261 SetAveragePt2(0.3); // GeV^2
>> 262 };
645 263
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 264
965 return; << 265 //G4cout<<"G4FTFParameters Out"<<G4endl;
966 } <<
967 266
>> 267 }
>> 268 //**********************************************************************************************
968 269