| 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 101902 2016-12-07 08:40:04Z gcosmo $
>> 28 // GEANT4 tag $Name: $
27 // 29 //
28 30
29 #include <utility> 31 #include <utility>
30 32
31 #include "G4FTFParameters.hh" 33 #include "G4FTFParameters.hh"
32 34
33 #include "G4ios.hh" 35 #include "G4ios.hh"
34 #include "G4PhysicalConstants.hh" 36 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
36 38
37 #include "G4ParticleDefinition.hh" 39 #include "G4ParticleDefinition.hh"
>> 40
38 #include "G4Proton.hh" 41 #include "G4Proton.hh"
39 #include "G4Neutron.hh" 42 #include "G4Neutron.hh"
>> 43
40 #include "G4PionPlus.hh" 44 #include "G4PionPlus.hh"
41 #include "G4PionMinus.hh" 45 #include "G4PionMinus.hh"
42 #include "G4KaonPlus.hh" 46 #include "G4KaonPlus.hh"
43 #include "G4KaonMinus.hh" 47 #include "G4KaonMinus.hh"
44 48
45 #include "G4CrossSectionDataSetRegistry.hh" <<
46 #include "G4VComponentCrossSection.hh" <<
47 #include "G4ComponentGGHadronNucleusXsc.hh" <<
48 #include "G4LundStringFragmentation.hh" <<
49 <<
50 #include "G4Exp.hh" 49 #include "G4Exp.hh"
51 #include "G4Log.hh" 50 #include "G4Log.hh"
52 #include "G4Pow.hh" 51 #include "G4Pow.hh"
53 52
54 #include "G4HadronicDeveloperParameters.hh" <<
55 #include "G4HadronicParameters.hh" <<
56 <<
57 //============================================ 53 //============================================================================
58 54
59 //#define debugFTFparams 55 //#define debugFTFparams
60 56
>> 57
61 //============================================ 58 //============================================================================
62 59
63 G4FTFParameters::G4FTFParameters() << 60 G4FTFParameters::G4FTFParameters() :
>> 61 FTFhNcmsEnergy( 0.0 ),
>> 62 FTFxsManager( 0 ),
>> 63 FTFXtotal( 0.0 ), FTFXelastic( 0.0 ), FTFXinelastic( 0.0 ), FTFXannihilation( 0.0 ),
>> 64 ProbabilityOfAnnihilation( 0.0 ), ProbabilityOfElasticScatt( 0.0 ),
>> 65 RadiusOfHNinteractions2( 0.0 ), FTFSlope( 0.0 ),
>> 66 AvaragePt2ofElasticScattering( 0.0 ), FTFGamma0( 0.0 ),
>> 67 DeltaProbAtQuarkExchange( 0.0 ), ProbOfSameQuarkExchange( 0.0 ),
>> 68 ProjMinDiffMass( 0.0 ), ProjMinNonDiffMass( 0.0 ), ProbLogDistrPrD(0.0),
>> 69 TarMinDiffMass( 0.0 ), TarMinNonDiffMass( 0.0 ),
>> 70 AveragePt2( 0.0 ), ProbLogDistr( 0.0 ),
>> 71 Pt2kink( 0.0 ),
>> 72 MaxNumberOfCollisions( 0.0 ), ProbOfInelInteraction( 0.0 ),
>> 73 CofNuclearDestructionPr( 0.0 ), CofNuclearDestruction( 0.0 ),
>> 74 R2ofNuclearDestruction( 0.0 ), ExcitationEnergyPerWoundedNucleon( 0.0 ),
>> 75 DofNuclearDestruction( 0.0 ), Pt2ofNuclearDestruction( 0.0 ), MaxPt2ofNuclearDestruction( 0.0 )
64 { 76 {
65 // Set-up alternative sets of FTF parameters << 77 for ( G4int i = 0; i < 4; i++ ) {
66 // Note that the very first tune (with index << 78 for ( G4int j = 0; j < 7; j++ ) {
67 // set of parameters, which does not need to << 79 ProcParams[i][j] = 0.0;
68 // the for loop below starts from 1 and not << 80 }
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 } 81 }
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 } 82 }
93 83
>> 84
>> 85 //============================================================================
>> 86
>> 87 G4FTFParameters::~G4FTFParameters() {}
>> 88
>> 89
94 //============================================ 90 //============================================================================
95 91
96 void G4FTFParameters::InitForInteraction( cons << 92 G4ThreadLocal bool G4FTFParameters::chipsComponentXSisInitialized = false;
97 G4in << 93 G4ThreadLocal G4ChipsComponentXS* G4FTFParameters::chipsComponentXSinstance = 0;
>> 94
>> 95
>> 96 //============================================================================
>> 97
>> 98 G4FTFParameters::G4FTFParameters( const G4ParticleDefinition* particle,
>> 99 G4int theA, G4int theZ, G4double PlabPerParticle ) :
>> 100 FTFhNcmsEnergy( 0.0 ),
>> 101 FTFxsManager( 0 ),
>> 102 FTFXtotal( 0.0 ), FTFXelastic( 0.0 ), FTFXinelastic( 0.0 ), FTFXannihilation( 0.0 ),
>> 103 ProbabilityOfAnnihilation( 0.0 ), ProbabilityOfElasticScatt( 0.0 ),
>> 104 RadiusOfHNinteractions2( 0.0 ), FTFSlope( 0.0 ),
>> 105 AvaragePt2ofElasticScattering( 0.0 ), FTFGamma0( 0.0 ),
>> 106 DeltaProbAtQuarkExchange( 0.0 ), ProbOfSameQuarkExchange( 0.0 ),
>> 107 ProjMinDiffMass( 0.0 ), ProjMinNonDiffMass( 0.0 ), ProbLogDistrPrD(0.0),
>> 108 TarMinDiffMass( 0.0 ), TarMinNonDiffMass( 0.0 ),
>> 109 AveragePt2( 0.0 ), ProbLogDistr( 0.0 ),
>> 110 Pt2kink( 0.0 ),
>> 111 MaxNumberOfCollisions( 0.0 ), ProbOfInelInteraction( 0.0 ),
>> 112 CofNuclearDestructionPr( 0.0 ), CofNuclearDestruction( 0.0 ),
>> 113 R2ofNuclearDestruction( 0.0 ), ExcitationEnergyPerWoundedNucleon( 0.0 ),
>> 114 DofNuclearDestruction( 0.0 ), Pt2ofNuclearDestruction( 0.0 ), MaxPt2ofNuclearDestruction( 0.0 )
98 { 115 {
99 Reset(); << 116 for ( G4int i = 0; i < 4; i++ ) {
>> 117 for ( G4int j = 0; j < 7; j++ ) {
>> 118 ProcParams[i][j] = 0.0;
>> 119 }
>> 120 }
100 121
101 G4int ProjectilePDGcode = particle->Ge 122 G4int ProjectilePDGcode = particle->GetPDGEncoding();
102 G4int ProjectileabsPDGcode = std::abs( Pr 123 G4int ProjectileabsPDGcode = std::abs( ProjectilePDGcode );
103 G4double ProjectileMass = particle->Ge 124 G4double ProjectileMass = particle->GetPDGMass();
104 G4double ProjectileMass2 = ProjectileMa 125 G4double ProjectileMass2 = ProjectileMass * ProjectileMass;
105 126
106 G4int ProjectileBaryonNumber( 0 ), AbsProjec 127 G4int ProjectileBaryonNumber( 0 ), AbsProjectileBaryonNumber( 0 ), AbsProjectileCharge( 0 );
107 G4bool ProjectileIsNucleus = false; 128 G4bool ProjectileIsNucleus = false;
108 << 129
109 if ( std::abs( particle->GetBaryonNumber() ) 130 if ( std::abs( particle->GetBaryonNumber() ) > 1 ) { // The projectile is a nucleus
110 ProjectileIsNucleus = true; 131 ProjectileIsNucleus = true;
111 ProjectileBaryonNumber = particle->GetB 132 ProjectileBaryonNumber = particle->GetBaryonNumber();
112 AbsProjectileBaryonNumber = std::abs( Proj 133 AbsProjectileBaryonNumber = std::abs( ProjectileBaryonNumber );
113 AbsProjectileCharge = std::abs( G4in << 134 AbsProjectileCharge = G4int( particle->GetPDGCharge() );
114 if ( ProjectileBaryonNumber > 1 ) { 135 if ( ProjectileBaryonNumber > 1 ) {
115 ProjectilePDGcode = 2212; ProjectileabsP 136 ProjectilePDGcode = 2212; ProjectileabsPDGcode = 2212; // Proton
116 } else { 137 } else {
117 ProjectilePDGcode = -2212; Projectileabs 138 ProjectilePDGcode = -2212; ProjectileabsPDGcode = 2212; // Anti-Proton
118 } 139 }
119 ProjectileMass = G4Proton::Proton()->GetP 140 ProjectileMass = G4Proton::Proton()->GetPDGMass();
120 ProjectileMass2 = sqr( ProjectileMass ); 141 ProjectileMass2 = sqr( ProjectileMass );
121 } 142 }
122 143
123 G4double TargetMass = G4Proton::Proton()->G 144 G4double TargetMass = G4Proton::Proton()->GetPDGMass();
124 G4double TargetMass2 = TargetMass * TargetMa 145 G4double TargetMass2 = TargetMass * TargetMass;
125 146
126 G4double Plab = PlabPerParticle; 147 G4double Plab = PlabPerParticle;
127 G4double Elab = std::sqrt( Plab*Plab + Proje 148 G4double Elab = std::sqrt( Plab*Plab + ProjectileMass2 );
128 G4double KineticEnergy = Elab - ProjectileMa 149 G4double KineticEnergy = Elab - ProjectileMass;
129 150
130 G4double S = ProjectileMass2 + TargetMass2 + 151 G4double S = ProjectileMass2 + TargetMass2 + 2.0*TargetMass*Elab;
131 152
132 #ifdef debugFTFparams 153 #ifdef debugFTFparams
133 G4cout << "--------- FTF Parameters -------- 154 G4cout << "--------- FTF Parameters --------------" << G4endl << "Proj Plab "
134 << ProjectilePDGcode << " " << Plab < 155 << ProjectilePDGcode << " " << Plab << G4endl << "Mass KinE " << ProjectileMass
135 << " " << KineticEnergy << G4endl << 156 << " " << KineticEnergy << G4endl << " A Z " << theA << " " << theZ << G4endl;
136 #endif 157 #endif
137 << 158
138 G4double Ylab, Xtotal( 0.0 ), Xelastic( 0.0 << 159 G4double Ylab, Xtotal, Xelastic, Xannihilation;
139 G4int NumberOfTargetNucleons; 160 G4int NumberOfTargetNucleons;
140 161
141 Ylab = 0.5 * G4Log( (Elab + Plab)/(Elab - Pl 162 Ylab = 0.5 * G4Log( (Elab + Plab)/(Elab - Plab) );
142 163
143 G4double ECMSsqr = S/GeV/GeV; 164 G4double ECMSsqr = S/GeV/GeV;
144 G4double SqrtS = std::sqrt( S )/GeV; 165 G4double SqrtS = std::sqrt( S )/GeV;
145 166
146 #ifdef debugFTFparams 167 #ifdef debugFTFparams
147 G4cout << "Sqrt(s) " << SqrtS << G4endl; 168 G4cout << "Sqrt(s) " << SqrtS << G4endl;
148 #endif 169 #endif
149 170
150 TargetMass /= GeV; TargetMass2 /= (G 171 TargetMass /= GeV; TargetMass2 /= (GeV*GeV);
151 ProjectileMass /= GeV; ProjectileMass2 /= (G 172 ProjectileMass /= GeV; ProjectileMass2 /= (GeV*GeV);
152 173
>> 174 // Andrea Dotti (13Jan2013):
>> 175 // The following lines are changed for G4MT. Originally the code was:
>> 176 // static G4ChipsComponentXS* _instance = new G4ChipsComponentXS(); // Witek Pokorski
>> 177 // Note the code could go back at original if _instance could be shared among threads
>> 178 if ( ! chipsComponentXSisInitialized ) {
>> 179 chipsComponentXSisInitialized = true;
>> 180 chipsComponentXSinstance = new G4ChipsComponentXS();
>> 181 }
>> 182 G4ChipsComponentXS* _instance = chipsComponentXSinstance;
>> 183 FTFxsManager = _instance;
>> 184
153 Plab /= GeV; 185 Plab /= GeV;
154 G4double Xftf = 0.0; 186 G4double Xftf = 0.0;
155 187
156 G4int NumberOfTargetProtons = theZ; 188 G4int NumberOfTargetProtons = theZ;
157 G4int NumberOfTargetNeutrons = theA - theZ; 189 G4int NumberOfTargetNeutrons = theA - theZ;
158 NumberOfTargetNucleons = NumberOfTargetProto 190 NumberOfTargetNucleons = NumberOfTargetProtons + NumberOfTargetNeutrons;
159 191
160 // ---------- hadron projectile ------------ << 192 if ( ProjectilePDGcode == 2212 || ProjectilePDGcode == 2112 ) { // Projectile is nucleon
161 if ( AbsProjectileBaryonNumber <= 1 ) { // << 193 G4ParticleDefinition* Proton = G4Proton::Proton(); //ALB
162 << 194 G4double XtotPP = FTFxsManager->GetTotalElementCrossSection( Proton, KineticEnergy, 1, 0 ); //ALB
163 // Interaction on P << 195
164 G4double xTtP = csGGinstance->GetTotalIsot << 196 G4ParticleDefinition* Neutron = G4Neutron::Neutron();
165 G4double xElP = csGGinstance->GetElasticIs << 197 G4double XtotPN = FTFxsManager->GetTotalElementCrossSection( Neutron, KineticEnergy, 1, 0 ); //ALB
166 << 198 G4double XelPP = FTFxsManager->GetElasticElementCrossSection( Proton, KineticEnergy, 1, 0 );
167 // Interaction on N << 199 G4double XelPN = FTFxsManager->GetElasticElementCrossSection( Neutron, KineticEnergy, 1, 0 );
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 200
179 #ifdef debugFTFparams 201 #ifdef debugFTFparams
180 G4cout<<"Estimated cross sections (total a << 202 G4cout << "XsPP " << XtotPP/millibarn << " " << XelPP/millibarn << G4endl
>> 203 << "XsPN " << XtotPN/millibarn << " " << XelPN/millibarn << G4endl;
181 #endif 204 #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 205
201 #ifdef debugFTFparams << 206 if ( ! ProjectileIsNucleus ) { // Projectile is hadron
202 G4cout << "XsPP (total and elastic) " << X << 207 Xtotal = ( NumberOfTargetProtons * XtotPP + NumberOfTargetNeutrons * XtotPN ) /
203 << "XsPN (total and elastic) " << X << 208 NumberOfTargetNucleons;
204 #endif << 209 Xelastic = ( NumberOfTargetProtons * XelPP + NumberOfTargetNeutrons * XelPN ) /
205 << 210 NumberOfTargetNucleons;
206 Xtotal = ( << 211 } else { // Projectile is a nucleus
>> 212 Xtotal = (
207 AbsProjectileCharge * Number 213 AbsProjectileCharge * NumberOfTargetProtons * XtotPP +
208 ( AbsProjectileBaryonNumber 214 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
209 NumberOfTargetNeutrons * 215 NumberOfTargetNeutrons * XtotPP
210 + 216 +
211 ( AbsProjectileCharge * Numb 217 ( AbsProjectileCharge * NumberOfTargetNeutrons +
212 ( AbsProjectileBaryonNumbe 218 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
213 NumberOfTargetProtons 219 NumberOfTargetProtons ) * XtotPN
214 ) / ( AbsProjectileBaryonNumbe 220 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
215 Xelastic= ( << 221 Xelastic= (
216 AbsProjectileCharge * Number 222 AbsProjectileCharge * NumberOfTargetProtons * XelPP +
217 ( AbsProjectileBaryonNumber 223 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
218 NumberOfTargetNeutrons * 224 NumberOfTargetNeutrons * XelPP
219 + 225 +
220 ( AbsProjectileCharge * Numb 226 ( AbsProjectileCharge * NumberOfTargetNeutrons +
221 ( AbsProjectileBaryonNumbe 227 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
222 NumberOfTargetProtons 228 NumberOfTargetProtons ) * XelPN
223 ) / ( AbsProjectileBaryonNumbe 229 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
>> 230 }
224 231
225 Xannihilation = 0.0; 232 Xannihilation = 0.0;
226 Xtotal /= millibarn; 233 Xtotal /= millibarn;
227 Xelastic /= millibarn; 234 Xelastic /= millibarn;
228 } <<
229 235
230 // ---------- The projectile is anti-baryon << 236 } else if ( ProjectilePDGcode < -1000 ) { // 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 237
240 G4double X_a( 0.0 ), X_b( 0.0 ), X_c( 0.0 238 G4double X_a( 0.0 ), X_b( 0.0 ), X_c( 0.0 ), X_d( 0.0 );
241 G4double MesonProdThreshold = ProjectileMa 239 G4double MesonProdThreshold = ProjectileMass + TargetMass +
242 ( 2.0 * 0.14 240 ( 2.0 * 0.14 + 0.016 ); // 2 Mpi + DeltaE;
243 241
244 if ( PlabPerParticle < 40.0*MeV ) { // Low 242 if ( PlabPerParticle < 40.0*MeV ) { // Low energy limits. Projectile at rest.
245 Xtotal = 1512.9; // mb 243 Xtotal = 1512.9; // mb
246 Xelastic = 473.2; // mb 244 Xelastic = 473.2; // mb
247 X_a = 625.1; // mb 245 X_a = 625.1; // mb
248 X_b = 9.780; // mb 246 X_b = 9.780; // mb
249 X_c = 49.989; // mb 247 X_c = 49.989; // mb
250 X_d = 6.614; // mb 248 X_d = 6.614; // mb
251 } else { // Total and elastic cross sectio 249 } else { // Total and elastic cross section of PbarP interactions a'la Arkhipov
252 G4double LogS = G4Log( ECMSsqr / 33.0625 250 G4double LogS = G4Log( ECMSsqr / 33.0625 );
253 G4double Xasmpt = 36.04 + 0.304*LogS*Log 251 G4double Xasmpt = 36.04 + 0.304*LogS*LogS; // mb
254 LogS = G4Log( SqrtS / 20.74 ); 252 LogS = G4Log( SqrtS / 20.74 );
255 G4double Basmpt = 11.92 + 0.3036*LogS*Lo 253 G4double Basmpt = 11.92 + 0.3036*LogS*LogS; // GeV^(-2)
256 G4double R0 = std::sqrt( 0.40874044*Xasm 254 G4double R0 = std::sqrt( 0.40874044*Xasmpt - Basmpt ); // GeV^(-1)
257 255
258 G4double FlowF = SqrtS / std::sqrt( ECMS 256 G4double FlowF = SqrtS / std::sqrt( ECMSsqr*ECMSsqr + ProjectileMass2*ProjectileMass2 +
259 Targ 257 TargetMass2*TargetMass2 - 2.0*ECMSsqr*ProjectileMass2
260 - 2. 258 - 2.0*ECMSsqr*TargetMass2
261 - 2. 259 - 2.0*ProjectileMass2*TargetMass2 );
262 260
263 Xtotal = Xasmpt * ( 1.0 + 13.55*FlowF/R0 261 Xtotal = Xasmpt * ( 1.0 + 13.55*FlowF/R0/R0/R0*
264 (1.0 - 4.47/Sq 262 (1.0 - 4.47/SqrtS + 12.38/ECMSsqr - 12.43/SqrtS/ECMSsqr) ); // mb
265 263
266 Xasmpt = 4.4 + 0.101*LogS*LogS; // mb 264 Xasmpt = 4.4 + 0.101*LogS*LogS; // mb
267 Xelastic = Xasmpt * ( 1.0 + 59.27*FlowF/ 265 Xelastic = Xasmpt * ( 1.0 + 59.27*FlowF/R0/R0/R0*
268 (1.0 - 6.95/ 266 (1.0 - 6.95/SqrtS + 23.54/ECMSsqr - 25.34/SqrtS/ECMSsqr ) ); // mb
269 267
>> 268 //G4cout << "Param Xtotal Xelastic " << Xtotal << " " << Xelastic << G4endl
>> 269 // << "FlowF " << FlowF << " SqrtS " << SqrtS << G4endl
>> 270 // << "Param Xelastic-NaN " << Xelastic << " "
>> 271 // << 1.5*16.654/pow(ECMSsqr/2.176/2.176,2.2) << " " << ECMSsqr << G4endl;
>> 272
270 X_a = 25.0*FlowF; // mb, 3-shirts diagr 273 X_a = 25.0*FlowF; // mb, 3-shirts diagram
271 274
272 if ( SqrtS < MesonProdThreshold ) { 275 if ( SqrtS < MesonProdThreshold ) {
273 X_b = 3.13 + 140.0*G4Pow::GetInstance( 276 X_b = 3.13 + 140.0*G4Pow::GetInstance()->powA( MesonProdThreshold - SqrtS, 2.5 ); // mb anti-quark-quark annihilation
274 Xelastic -= 3.0*X_b; // Xel-X(PbarP-> 277 Xelastic -= 3.0*X_b; // Xel-X(PbarP->NNbar)
275 } else { 278 } else {
276 X_b = 6.8/SqrtS; // mb anti-quark-qua 279 X_b = 6.8/SqrtS; // mb anti-quark-quark annihilation
277 Xelastic -= 3.0*X_b; // Xel-X(PbarP-> 280 Xelastic -= 3.0*X_b; // Xel-X(PbarP->NNbar)
278 } 281 }
279 282
280 X_c = 2.0*FlowF*sqr( ProjectileMass + Ta 283 X_c = 2.0*FlowF*sqr( ProjectileMass + TargetMass )/ECMSsqr; // mb rearrangement
281 284
282 X_d = 23.3/ECMSsqr; // mb anti-quark-qu 285 X_d = 23.3/ECMSsqr; // mb anti-quark-quark string creation
283 } 286 }
284 287
>> 288 //G4cout << "Param Xtotal Xelastic " << Xtotal << " " << Xelastic << G4endl
>> 289 // << "Para a b c d " << X_a << " " << X_b << " " << X_c << " " << X_d << G4endl;
>> 290 // << "Para a b c d " << X_a << " " << 5.*X_b << " " << 5.*X_c << " " << 6.*X_d
>> 291 // << G4endl;
>> 292
285 G4double Xann_on_P( 0.0), Xann_on_N( 0.0 ) 293 G4double Xann_on_P( 0.0), Xann_on_N( 0.0 );
286 294
287 if ( ProjectilePDGcode == -2212 ) { << 295 if ( ProjectilePDGcode == -2212 ) { // Pbar+P/N
288 Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_ 296 Xann_on_P = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0;
289 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_ 297 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
290 } else if ( ProjectilePDGcode == -2112 ) { << 298 } else if ( ProjectilePDGcode == -2112 ) { // NeutrBar+P/N
291 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_ 299 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*4.0;
292 Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_ 300 Xann_on_N = X_a + X_b*5.0 + X_c*5.0 + X_d*6.0;
293 } else if ( ProjectilePDGcode == -3122 ) { << 301 } else if ( ProjectilePDGcode == -3122 ) { // LambdaBar+P/N
294 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_ 302 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
295 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_ 303 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
296 } else if ( ProjectilePDGcode == -3112 ) { << 304 } else if ( ProjectilePDGcode == -3112 ) { // Sigma-Bar+P/N
297 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_ 305 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
298 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_ 306 Xann_on_N = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
299 } else if ( ProjectilePDGcode == -3212 ) { << 307 } else if ( ProjectilePDGcode == -3212 ) { // Sigma0Bar+P/N
300 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_ 308 Xann_on_P = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
301 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_ 309 Xann_on_N = X_a + X_b*3.0 + X_c*3.0 + X_d*2.0;
302 } else if ( ProjectilePDGcode == -3222 ) { << 310 } else if ( ProjectilePDGcode == -3222 ) { // Sigma+Bar+P/N
303 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_ 311 Xann_on_P = X_a + X_b*4.0 + X_c*4.0 + X_d*2.0;
304 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_ 312 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
305 } else if ( ProjectilePDGcode == -3312 ) { << 313 } else if ( ProjectilePDGcode == -3312 ) { // Xi-Bar+P/N
306 Xann_on_P = X_a + X_b*1.0 + X_c*1.0 + X_ 314 Xann_on_P = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
307 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_ 315 Xann_on_N = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
308 } else if ( ProjectilePDGcode == -3322 ) { << 316 } else if ( ProjectilePDGcode == -3322 ) { // Xi0Bar+P/N
309 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_ 317 Xann_on_P = X_a + X_b*2.0 + X_c*2.0 + X_d*0.0;
310 Xann_on_N = X_a + X_b*1.0 + X_c*1.0 + X_ 318 Xann_on_N = X_a + X_b*1.0 + X_c*1.0 + X_d*0.0;
311 } else if ( ProjectilePDGcode == -3334 ) { << 319 } else if ( ProjectilePDGcode == -3334 ) { // Omega-Bar+P/N
312 Xann_on_P = X_a + X_b*0.0 + X_c*0.0 + X_ 320 Xann_on_P = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
313 Xann_on_N = X_a + X_b*0.0 + X_c*0.0 + X_ 321 Xann_on_N = X_a + X_b*0.0 + X_c*0.0 + X_d*0.0;
314 } else { 322 } else {
315 G4cout << "Unknown anti-baryon for FTF a 323 G4cout << "Unknown anti-baryon for FTF annihilation" << G4endl;
316 } 324 }
317 325
318 //G4cout << "Sum " << Xann_on_P < 326 //G4cout << "Sum " << Xann_on_P << G4endl;
319 327
320 if ( ! ProjectileIsNucleus ) { // Project 328 if ( ! ProjectileIsNucleus ) { // Projectile is anti-baryon
321 Xannihilation = ( NumberOfTargetProtons 329 Xannihilation = ( NumberOfTargetProtons * Xann_on_P + NumberOfTargetNeutrons * Xann_on_N )
322 / NumberOfTargetNucleons 330 / NumberOfTargetNucleons;
323 } else { // Projectile is a nucleus 331 } else { // Projectile is a nucleus
324 Xannihilation = ( 332 Xannihilation = (
325 ( AbsProjectileCharge 333 ( AbsProjectileCharge * NumberOfTargetProtons +
326 ( AbsProjectileBaryo 334 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
327 NumberOfTargetNeutro 335 NumberOfTargetNeutrons ) * Xann_on_P
328 + 336 +
329 ( AbsProjectileCharge 337 ( AbsProjectileCharge * NumberOfTargetNeutrons +
330 ( AbsProjectileBaryo 338 ( AbsProjectileBaryonNumber - AbsProjectileCharge ) *
331 NumberOfTargetProton 339 NumberOfTargetProtons ) * Xann_on_N
332 ) / ( AbsProjectileBaryo 340 ) / ( AbsProjectileBaryonNumber * NumberOfTargetNucleons );
333 } 341 }
334 342
335 //G4double Xftf = 0.0; 343 //G4double Xftf = 0.0;
336 MesonProdThreshold = ProjectileMass + Targ 344 MesonProdThreshold = ProjectileMass + TargetMass + (0.14 + 0.08); // Mpi + DeltaE
337 if ( SqrtS > MesonProdThreshold ) { 345 if ( SqrtS > MesonProdThreshold ) {
338 Xftf = 36.0 * ( 1.0 - MesonProdThreshold 346 Xftf = 36.0 * ( 1.0 - MesonProdThreshold/SqrtS );
339 } 347 }
340 348
341 Xtotal = Xelastic + Xannihilation + Xftf; 349 Xtotal = Xelastic + Xannihilation + Xftf;
342 350
343 #ifdef debugFTFparams 351 #ifdef debugFTFparams
344 G4cout << "Plab Xtotal, Xelastic Xinel Xf 352 G4cout << "Plab Xtotal, Xelastic Xinel Xftf " << Plab << " " << Xtotal << " " << Xelastic
345 << " " << Xtotal - Xelastic << " " << 353 << " " << Xtotal - Xelastic << " " << Xtotal - Xelastic - Xannihilation << G4endl
346 << "Plab Xelastic/Xtotal, Xann/Xin 354 << "Plab Xelastic/Xtotal, Xann/Xin " << Plab << " " << Xelastic/Xtotal << " "
347 << Xannihilation/(Xtotal - Xelastic 355 << Xannihilation/(Xtotal - Xelastic) << G4endl;
348 #endif 356 #endif
349 357
350 } << 358 } else if ( ProjectilePDGcode == 211 ) { // Projectile is PionPlus
>> 359
>> 360 G4double XtotPiP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 361 G4ParticleDefinition* PionMinus = G4PionMinus::PionMinus();
>> 362 G4double XtotPiN = FTFxsManager->GetTotalElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
>> 363 G4double XelPiP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 364 G4double XelPiN = FTFxsManager->GetElasticElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
>> 365 Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
>> 366 / NumberOfTargetNucleons;
>> 367 Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
>> 368 / NumberOfTargetNucleons;
>> 369 Xannihilation = 0.0;
>> 370 Xtotal /= millibarn;
>> 371 Xelastic /= millibarn;
>> 372
>> 373 } else if ( ProjectilePDGcode == -211 ) { // Projectile is PionMinus
>> 374
>> 375 G4double XtotPiP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 376 G4ParticleDefinition* PionPlus = G4PionPlus::PionPlus();
>> 377 G4double XtotPiN = FTFxsManager->GetTotalElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
>> 378 G4double XelPiP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 379 G4double XelPiN = FTFxsManager->GetElasticElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
>> 380 Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
>> 381 / NumberOfTargetNucleons;
>> 382 Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
>> 383 / NumberOfTargetNucleons;
>> 384 Xannihilation = 0.0;
>> 385 Xtotal /= millibarn;
>> 386 Xelastic /= millibarn;
>> 387
>> 388 } else if ( ProjectilePDGcode == 111 ) { // Projectile is PionZero
>> 389
>> 390 G4ParticleDefinition* PionPlus = G4PionPlus::PionPlus();
>> 391 G4double XtotPipP = FTFxsManager->GetTotalElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
>> 392 G4ParticleDefinition* PionMinus = G4PionMinus::PionMinus();
>> 393 G4double XtotPimP = FTFxsManager->GetTotalElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
>> 394 G4double XelPipP = FTFxsManager->GetElasticElementCrossSection( PionPlus, KineticEnergy, 1, 0 );
>> 395 G4double XelPimP = FTFxsManager->GetElasticElementCrossSection( PionMinus, KineticEnergy, 1, 0 );
>> 396 G4double XtotPiP = ( XtotPipP + XtotPimP ) / 2.0;
>> 397 G4double XtotPiN = XtotPiP;
>> 398 G4double XelPiP = ( XelPipP + XelPimP ) / 2.0;
>> 399 G4double XelPiN = XelPiP;
>> 400 Xtotal = ( NumberOfTargetProtons * XtotPiP + NumberOfTargetNeutrons * XtotPiN )
>> 401 / NumberOfTargetNucleons;
>> 402 Xelastic = ( NumberOfTargetProtons * XelPiP + NumberOfTargetNeutrons * XelPiN )
>> 403 / NumberOfTargetNucleons;
>> 404 Xannihilation = 0.0;
>> 405 Xtotal /= millibarn;
>> 406 Xelastic /= millibarn;
>> 407
>> 408 } else if ( ProjectilePDGcode == 321 ) { // Projectile is KaonPlus
>> 409
>> 410 G4double XtotKP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 411 G4ParticleDefinition* KaonMinus = G4KaonMinus::KaonMinus();
>> 412 G4double XtotKN = FTFxsManager->GetTotalElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
>> 413 G4double XelKP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 414 G4double XelKN = FTFxsManager->GetElasticElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
>> 415 Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
>> 416 / NumberOfTargetNucleons;
>> 417 Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
>> 418 / NumberOfTargetNucleons;
>> 419 Xannihilation = 0.0;
>> 420 Xtotal /= millibarn;
>> 421 Xelastic /= millibarn;
>> 422
>> 423 } else if ( ProjectilePDGcode == -321 ) { // Projectile is KaonMinus
351 424
352 if ( Xtotal == 0.0 ) { // Projectile is und << 425 G4double XtotKP = FTFxsManager->GetTotalElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 426 G4ParticleDefinition* KaonPlus = G4KaonPlus::KaonPlus();
>> 427 G4double XtotKN = FTFxsManager->GetTotalElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
>> 428 G4double XelKP = FTFxsManager->GetElasticElementCrossSection( particle, KineticEnergy, 1, 0 );
>> 429 G4double XelKN = FTFxsManager->GetElasticElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
>> 430 Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
>> 431 / NumberOfTargetNucleons;
>> 432 Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
>> 433 / NumberOfTargetNucleons;
>> 434 Xannihilation = 0.0;
>> 435 Xtotal /= millibarn;
>> 436 Xelastic /= millibarn;
>> 437
>> 438 } else if ( ProjectilePDGcode == 311 || ProjectilePDGcode == 130 ||
>> 439 ProjectilePDGcode == 310 ) { // Projectile is KaonZero
>> 440
>> 441 G4ParticleDefinition* KaonPlus = G4KaonPlus::KaonPlus();
>> 442 G4double XtotKpP = FTFxsManager->GetTotalElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
>> 443 G4ParticleDefinition* KaonMinus = G4KaonMinus::KaonMinus();
>> 444 G4double XtotKmP = FTFxsManager->GetTotalElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
>> 445 G4double XelKpP = FTFxsManager->GetElasticElementCrossSection( KaonPlus, KineticEnergy, 1, 0 );
>> 446 G4double XelKmP = FTFxsManager->GetElasticElementCrossSection( KaonMinus, KineticEnergy, 1, 0 );
>> 447 G4double XtotKP = ( XtotKpP + XtotKmP ) / 2.0;
>> 448 G4double XtotKN = XtotKP;
>> 449 G4double XelKP = ( XelKpP + XelKmP ) / 2.0;
>> 450 G4double XelKN = XelKP;
>> 451 Xtotal = ( NumberOfTargetProtons * XtotKP + NumberOfTargetNeutrons * XtotKN )
>> 452 / NumberOfTargetNucleons;
>> 453 Xelastic = ( NumberOfTargetProtons * XelKP + NumberOfTargetNeutrons * XelKN )
>> 454 / NumberOfTargetNucleons;
>> 455 Xannihilation = 0.0;
>> 456 Xtotal /= millibarn;
>> 457 Xelastic /= millibarn;
353 458
354 const G4ParticleDefinition* Proton = G4Pro << 459 } else { // Projectile is undefined, Nucleon assumed
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 460
>> 461 G4ParticleDefinition* Proton = G4Proton::Proton();
>> 462 G4double XtotPP = FTFxsManager->GetTotalElementCrossSection( Proton, KineticEnergy, 1, 0 );
>> 463 G4ParticleDefinition* Neutron = G4Neutron::Neutron();
>> 464 G4double XtotPN = FTFxsManager->GetTotalElementCrossSection( Neutron, KineticEnergy, 1, 0 );
>> 465 G4double XelPP = FTFxsManager->GetElasticElementCrossSection( Proton, KineticEnergy, 1, 0 );
>> 466 G4double XelPN = FTFxsManager->GetElasticElementCrossSection( Neutron, KineticEnergy, 1, 0 );
363 Xtotal = ( NumberOfTargetProtons * Xtot 467 Xtotal = ( NumberOfTargetProtons * XtotPP + NumberOfTargetNeutrons * XtotPN )
364 / NumberOfTargetNucleons; 468 / NumberOfTargetNucleons;
365 Xelastic = ( NumberOfTargetProtons * XelP 469 Xelastic = ( NumberOfTargetProtons * XelPP + NumberOfTargetNeutrons * XelPN )
366 / NumberOfTargetNucleons; 470 / NumberOfTargetNucleons;
367 Xannihilation = 0.0; 471 Xannihilation = 0.0;
368 Xtotal /= millibarn; 472 Xtotal /= millibarn;
369 Xelastic /= millibarn; 473 Xelastic /= millibarn;
>> 474
370 }; 475 };
371 476
372 // Geometrical parameters 477 // Geometrical parameters
373 SetTotalCrossSection( Xtotal ); 478 SetTotalCrossSection( Xtotal );
374 SetElasticCrossSection( Xelastic ); << 479 SetElastisCrossSection( Xelastic );
375 SetInelasticCrossSection( Xtotal - Xelastic 480 SetInelasticCrossSection( Xtotal - Xelastic );
376 481
>> 482 //G4cout << "Plab Xtotal, Xelastic Xinel Xftf " << Plab << " " << Xtotal << " " << Xelastic
>> 483 // << " " << Xtotal - Xelastic << " " << Xtotal - Xelastic - Xannihilation << G4endl;
>> 484 //if (Xtotal - Xelastic != 0.0 ) {
>> 485 // G4cout << "Plab Xelastic/Xtotal, Xann/Xin " << Plab << " " << Xelastic/Xtotal
>> 486 // << " " << Xannihilation / (Xtotal - Xelastic) << G4endl;
>> 487 //} else {
>> 488 // G4cout << "Plab Xelastic/Xtotal, Xann " << Plab << " " << Xelastic/Xtotal
>> 489 // << " " << Xannihilation << G4endl;
>> 490 //}
>> 491 //G4int Uzhi; G4cin >> Uzhi;
>> 492
377 // Interactions with elastic and inelastic c 493 // Interactions with elastic and inelastic collisions
378 SetProbabilityOfElasticScatt( Xtotal, Xelast 494 SetProbabilityOfElasticScatt( Xtotal, Xelastic );
379 <<
380 SetRadiusOfHNinteractions2( Xtotal/pi/10.0 ) 495 SetRadiusOfHNinteractions2( Xtotal/pi/10.0 );
381 << 496 if ( Xtotal - Xelastic == 0.0 ) {
382 if ( ( Xtotal - Xelastic ) == 0.0 ) { <<
383 SetProbabilityOfAnnihilation( 0.0 ); 497 SetProbabilityOfAnnihilation( 0.0 );
384 } else { 498 } else {
385 SetProbabilityOfAnnihilation( Xannihilatio 499 SetProbabilityOfAnnihilation( Xannihilation / (Xtotal - Xelastic) );
386 } 500 }
387 501
388 if(Xelastic > 0.0) { << 502 // No elastic scattering
389 SetSlope( Xtotal*Xtotal/16.0/pi/Xelastic/0 << 503 //SetProbabilityOfElasticScatt( Xtotal, 0.0 );
390 << 504 //SetRadiusOfHNinteractions2( (Xtotal - Xelastic)/pi/10.0 );
391 // Parameters of elastic scattering << 505 //SetProbabilityOfAnnihilation( 1.0 );
392 // Gaussian parametrization of elastic sca << 506 //SetProbabilityOfAnnihilation( 0.0 );
393 SetAvaragePt2ofElasticScattering( 1.0/( Xt << 507
394 } else { << 508 SetSlope( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 ); // Slope parameter of elastic scattering
395 SetSlope(1.0); << 509 // (GeV/c)^(-2))
396 SetAvaragePt2ofElasticScattering( 0.0); << 510 //G4cout << "Slope " << GetSlope() << G4endl;
397 } <<
398 SetGamma0( GetSlope()*Xtotal/10.0/2.0/pi ); 511 SetGamma0( GetSlope()*Xtotal/10.0/2.0/pi );
399 512
400 G4double Xinel = Xtotal - Xelastic; << 513 // Parameters of elastic scattering
>> 514 // Gaussian parametrization of elastic scattering amplitude assumed
>> 515 SetAvaragePt2ofElasticScattering( 1.0/( Xtotal*Xtotal/16.0/pi/Xelastic/0.3894 )*GeV*GeV );
>> 516 //G4cout << "AvaragePt2ofElasticScattering " << GetAvaragePt2ofElasticScattering() << G4endl;
401 517
402 #ifdef debugFTFparams << 518 // Parameters of excitations
403 G4cout<< "Slope of hN elastic scattering" << <<
404 G4cout << "AvaragePt2ofElasticScattering " < <<
405 G4cout<<"Parameters of excitation for projec <<
406 #endif <<
407 519
408 if ( (ProjectilePDGcode == 2212) || (Project << 520 G4double Xinel = Xtotal - Xelastic;
409 << 521 //G4cout << "Param ProjectilePDGcode " << ProjectilePDGcode << G4endl;
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 522
>> 523 if ( ProjectilePDGcode > 1000 ) { // Projectile is baryon
415 // Proc# A1 B1 A2 524 // Proc# A1 B1 A2 B2 A3 Atop Ymin
416 /* original hadr-string-diff-V10-03-07 (si <<
417 SetParams( 0, 13.71, 1.75, -2 525 SetParams( 0, 13.71, 1.75, -214.5, 4.25, 0.0, 0.5 , 1.1 ); // Qexchange without Exc.
418 SetParams( 1, 25.0, 1.0, -5 526 SetParams( 1, 25.0, 1.0, -50.34, 1.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc.
419 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*1 << 527 if( Xinel > 0.) {
420 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*1 << 528 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93);// Projectile diffraction
421 SetParams( 4, 1.0, 0.0 , -2 << 529 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93);// Target diffraction
422 */ << 530 SetParams( 4, 1.0, 0.0 , -2.01 , 0.5 , 0.0, 0.0 , 1.4 );// Qexchange with Exc. Additional multiply
423 << 531 } else {
424 // Proc# << 532 SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
425 SetParams( 0, fArrayParCollBaryonProj[inde << 533 SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
426 fArrayParCollBaryonProj[indexTune] << 534 SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
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 } 535 }
455 536
456 if ( (AbsProjectileBaryonNumber > 10 || Nu << 537 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
457 // It is not decided what to do with dif 538 // It is not decided what to do with diffraction dissociation in Had-Nucl and Nucl-Nucl interactions
458 // For the moment both ProjDiffDisso & T << 539 SetParams( 2, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
459 // so both projectile and target diffrac << 540 //SetParams( 3, 0.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
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 } 541 }
465 542
466 SetDeltaProbAtQuarkExchange( fArrayParColl << 543 SetDeltaProbAtQuarkExchange( 0.0 );
467 <<
468 if ( NumberOfTargetNucleons > 26 ) { 544 if ( NumberOfTargetNucleons > 26 ) {
469 SetProbOfSameQuarkExchange( 1.0 ); << 545 SetProbOfSameQuarkExchange( 1.0);
470 } else { 546 } else {
471 SetProbOfSameQuarkExchange( fArrayParCol << 547 SetProbOfSameQuarkExchange( 0.0 );
472 } 548 }
>> 549 SetProjMinDiffMass( 1.16 ); // GeV
>> 550 SetProjMinNonDiffMass( 1.16 ); // GeV
>> 551 SetTarMinDiffMass( 1.16 ); // GeV
>> 552 SetTarMinNonDiffMass( 1.16 ); // GeV
>> 553 SetAveragePt2( 0.15 ); // GeV^2
>> 554 SetProbLogDistrPrD( 0.3 ); // Before it was: 0.5
>> 555 SetProbLogDistr(0.3 ); // Before it was: 0.5
473 556
474 SetProjMinDiffMass( fArrayParCollBaryon << 557 } else if( ProjectilePDGcode < -1000 ) { // Projectile is anti_baryon
475 SetProjMinNonDiffMass( fArrayParCollBaryon <<
476 <<
477 SetTarMinDiffMass( fArrayParCollBaryon <<
478 SetTarMinNonDiffMass( fArrayParCollBaryon <<
479 558
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 559 // Proc# A1 B1 A2 B2 A3 Atop Ymin
493 SetParams( 0, 0.0 , 0.0 , 0 560 SetParams( 0, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 ); // Qexchange without Exc.
494 SetParams( 1, 0.0 , 0.0 , 0 561 SetParams( 1, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 1000.0 ); // Qexchange with Exc.
495 if ( Xinel > 0.) { << 562 if( Xinel > 0.) {
496 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel << 563 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Projectile diffraction
497 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel << 564 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Target diffraction
498 SetParams( 4, 1.0, 0.0 , << 565 SetParams( 4, 1.0, 0.0 , 0.0, 0.0 , 0.0, 0.0 , 0.93 ); // Qexchange with Exc. Additional multiply
499 } else { 566 } else {
500 SetParams( 2, 0.0, 0.0, 0.0, 0.0, 0.0, 0 << 567 SetParams( 2, 0.0, 0.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 << 568 SetParams( 3, 0.0, 0.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 << 569 SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
503 } 570 }
504 571
505 if ( AbsProjectileBaryonNumber > 10 || N << 572 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
506 // It is not decided what to do with dif << 573 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
507 // For the moment both ProjDiffDisso & T << 574 //SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
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 } 575 }
514 576
515 SetDeltaProbAtQuarkExchange( 0.0 ); 577 SetDeltaProbAtQuarkExchange( 0.0 );
516 SetProbOfSameQuarkExchange( 0.0 ); 578 SetProbOfSameQuarkExchange( 0.0 );
517 SetProjMinDiffMass( ProjectileMass + 0.22 579 SetProjMinDiffMass( ProjectileMass + 0.22 ); // GeV
518 SetProjMinNonDiffMass( ProjectileMass + 0. 580 SetProjMinNonDiffMass( ProjectileMass + 0.22 ); // GeV
519 SetTarMinDiffMass( TargetMass + 0.22 ); 581 SetTarMinDiffMass( TargetMass + 0.22 ); // GeV
520 SetTarMinNonDiffMass( TargetMass + 0.22 ); 582 SetTarMinNonDiffMass( TargetMass + 0.22 ); // GeV
521 SetAveragePt2( 0.3 ); << 583 SetAveragePt2( 0.15 ); // GeV^2
522 SetProbLogDistrPrD( 0.55 ); << 584 SetProbLogDistrPrD( 0.3 );
523 SetProbLogDistr( 0.55 ); << 585 SetProbLogDistr( 0.3 );
524 586
525 } else if ( ProjectileabsPDGcode == 211 || << 587 } else if ( ProjectileabsPDGcode == 211 || ProjectilePDGcode == 111 ) { // Projectile is Pion
526 << 588
527 const G4int indexTune = G4FTFTunings::Inst << 589 // Proc# A1 B1 A2 B2 A3 Atop Ymin
528 << 590 SetParams( 0, 720.0, 2.5 , 2.3 , 1.0, 0., 1. , 2.7 );
529 // Proc# A1 B1 A2 << 591 SetParams( 1, 12.87, 0.5 , -44.91, 1.0, 0., 0. , 2.5 );
530 /* --> original code << 592 SetParams( 2, 0.086, 0. , -0.3 , 0.5, 0., 0. , 2.5 );
531 SetParams( 0, 150.0, 1.8 , -247. << 593 SetParams( 3, 32.8, 1.0 , -114.5 , 1.5, 0.084, 0. , 2.5 );
532 SetParams( 1, 5.77, 0.6 , -5.7 << 594 SetParams( 4, 1.0, 0.0 , -3.49, 0.5, 0.0, 0. , 2.5 ); // Qexchange with Exc. Additional multiply
533 SetParams( 2, 2.27, 0.5 , -98052. << 595
534 SetParams( 3, 7.0, 0.9, -85.2 << 596 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
535 SetParams( 4, 1.0, 0.0 , -11.0 << 597 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
536 */ << 598 //SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
537 // Proc# << 599 }
538 SetParams( 0, fArrayParCollPionProj[indexT << 600
539 fArrayParCollPionProj[indexTune].G << 601 SetDeltaProbAtQuarkExchange( 0.56 ); // (0.35)
540 fArrayParCollPionProj[indexTune].GetProc << 602 SetProjMinDiffMass( 0.5 ); // (0.5) // GeV
541 fArrayParCollPionProj[indexTune].G << 603 SetProjMinNonDiffMass( 0.5 ); // (0.5) // GeV
542 fArrayParCollPionProj[indexTune].GetProc << 604 SetTarMinDiffMass( 1.16 ); // GeV
543 fArrayParCollPionProj[indexTune].G << 605 SetTarMinNonDiffMass( 1.16 ); // GeV
544 fArrayParCollPionProj[indexTune].GetProc << 606 SetAveragePt2( 0.15 ); // GeV^2
545 SetParams( 1, fArrayParCollPionProj[indexT << 607 SetProbLogDistrPrD( 0.3 );
546 fArrayParCollPionProj[indexTune].G << 608 SetProbLogDistr( 0.3 );
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 609
607 } else if ( ProjectileabsPDGcode == 321 || 610 } else if ( ProjectileabsPDGcode == 321 || ProjectileabsPDGcode == 311 ||
608 ProjectilePDGcode == 130 || 611 ProjectilePDGcode == 130 || ProjectilePDGcode == 310 ) { // Projectile is Kaon
609 612
610 // Proc# A1 B1 A2 613 // Proc# A1 B1 A2 B2 A3 Atop Ymin
611 SetParams( 0, 60.0 , 2.5 , 0 614 SetParams( 0, 60.0 , 2.5 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Qexchange without Exc.
612 SetParams( 1, 6.0 , 1.0 , -24 615 SetParams( 1, 6.0 , 1.0 , -24.33 , 2.0 , 0.0, 0.0 , 1.40 ); // Qexchange with Exc.
613 SetParams( 2, 2.76, 1.2 , -22 616 SetParams( 2, 2.76, 1.2 , -22.5 , 2.7 ,0.04, 0.0 , 1.40 ); // Projectile diffraction
614 SetParams( 3, 1.09, 0.5 , -8 617 SetParams( 3, 1.09, 0.5 , -8.88 , 2. ,0.05, 0.0 , 1.40 ); // Target diffraction
615 SetParams( 4, 1.0, 0.0 , 0 618 SetParams( 4, 1.0, 0.0 , 0.0 , 0.0 , 0.0, 0.0 , 0.93 ); // Qexchange with Exc. Additional multiply
616 if ( AbsProjectileBaryonNumber > 10 || N << 619
>> 620 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
617 SetParams( 2, 0.0 , 0.0 , 621 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
618 SetParams( 3, 0.0 , 0.0 , << 622 //SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
619 } 623 }
620 624
621 SetDeltaProbAtQuarkExchange( 0.6 ); 625 SetDeltaProbAtQuarkExchange( 0.6 );
622 SetProjMinDiffMass( 0.7 ); // GeV << 626 SetProjMinDiffMass( 0.7 ); // (1.4) // (0.7) // GeV
623 SetProjMinNonDiffMass( 0.7 ); // GeV << 627 SetProjMinNonDiffMass( 0.7 ); // (1.4) // (0.7) // GeV
624 SetTarMinDiffMass( 1.16 ); // GeV << 628 SetTarMinDiffMass( 1.16 ); // GeV
625 SetTarMinNonDiffMass( 1.16 ); // GeV << 629 SetTarMinNonDiffMass( 1.16 ); // GeV
626 SetAveragePt2( 0.3 ); // GeV^2 << 630 SetAveragePt2( 0.15 ); // GeV^2
627 SetProbLogDistrPrD( 0.55 ); << 631 SetProbLogDistrPrD( 0.5 );
628 SetProbLogDistr( 0.55 ); << 632 SetProbLogDistr( 0.3 );
629 <<
630 } else { // Projectile is not p, n, Pi0, Pi <<
631 <<
632 if ( ProjectileabsPDGcode > 1000 ) { // T <<
633 // Proc# A1 B1 <<
634 SetParams( 0, 13.71, 1.75, <<
635 SetParams( 1, 25.0, 1.0, - <<
636 if ( Xinel > 0.) { <<
637 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xin <<
638 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xin <<
639 SetParams( 4, 1.0, 0.0 , <<
640 } else { <<
641 SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0 <<
642 SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0 <<
643 SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0 <<
644 } <<
645 633
646 } else { // The projectile is a meson as << 634 } else { // Projectile is undefined, Nucleon assumed
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 635
655 if ( AbsProjectileBaryonNumber > 10 || N << 636 // Proc# A1 B1 A2 B2 A3 Atop Ymin
>> 637 SetParams( 0, 13.71, 1.75, -214.5, 4.25, 0.0, 0.5 , 1.1 ); // Qexchange without Exc.
>> 638 SetParams( 1, 25.0, 1.0, -50.34, 1.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc.
>> 639 if( Xinel > 0.) {
>> 640 SetParams( 2, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Projectile diffraction
>> 641 SetParams( 3, 6.0/Xinel, 0.0 ,-6.0/Xinel*16.28, 3.0 , 0.0, 0.0 , 0.93); // Target diffraction
>> 642 SetParams( 4, 1.0, 0.0 , -2.01 , 0.5 , 0.0, 0.0 , 1.4 ); // Qexchange with Exc. Additional multiply
>> 643 } else {
>> 644 SetParams( 2, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
>> 645 SetParams( 3, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
>> 646 SetParams( 4, 0.0, 0.0 ,0.0, 0.0 , 0.0, 0.0 , 0.0);
>> 647 }
>> 648 if ( AbsProjectileBaryonNumber > 1 || NumberOfTargetNucleons > 1 ) {
656 SetParams( 2, 0.0 , 0.0 , 649 SetParams( 2, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Projectile diffraction
657 SetParams( 3, 0.0 , 0.0 , << 650 //SetParams( 3, 0.0 , 0.0 , 0.0 , 0.0 , 0.0, 0.0 , -100.0 ); // Target diffraction
658 } 651 }
659 << 652 SetDeltaProbAtQuarkExchange( 0.0 ); // 7 June 2011
660 SetDeltaProbAtQuarkExchange( 0.0 ); <<
661 SetProbOfSameQuarkExchange( 0.0 ); 653 SetProbOfSameQuarkExchange( 0.0 );
662 << 654 SetProjMinDiffMass( ProjectileMass + 0.22 ); // GeV
663 SetProjMinDiffMass( GetMinMass(particle << 655 SetProjMinNonDiffMass( ProjectileMass + 0.22 ); // GeV
664 SetProjMinNonDiffMass( GetMinMass(particle << 656 SetTarMinDiffMass( TargetMass + 0.22 ); // GeV
665 << 657 SetTarMinNonDiffMass( TargetMass + 0.22 ); // GeV
666 const G4ParticleDefinition* Neutron = G4Ne << 658 SetAveragePt2( 0.15 ); // GeV^2
667 SetTarMinDiffMass( GetMinMass(Neutron)/ << 659 SetProbLogDistrPrD( 0.3 );
668 SetTarMinNonDiffMass( GetMinMass(Neutron)/ << 660 SetProbLogDistr( 0.3 );
669 <<
670 SetAveragePt2( 0.3 ); // GeV^2 <<
671 SetProbLogDistrPrD( 0.55 ); <<
672 SetProbLogDistr( 0.55 ); <<
673 661
674 } 662 }
675 663
676 #ifdef debugFTFparams << 664 // Set parameters of a string kink
677 G4cout<<"DeltaProbAtQuarkExchange "<< GetDel << 665 // SetPt2Kink( 6.0*GeV*GeV );
678 G4cout<<"ProbOfSameQuarkExchange "<< GetPro << 666 SetPt2Kink( 0.0*GeV*GeV ); // Uzhi Oct 2014 to switch off kinky strings
679 G4cout<<"ProjMinDiffMass "<< GetPro << 667 G4double Puubar( 1.0/3.0 ), Pddbar( 1.0/3.0 ), Pssbar( 1.0/3.0 ); // SU(3) symmetry
680 G4cout<<"ProjMinNonDiffMass "<< GetPro << 668 //G4double Puubar( 0.41 ), Pddbar( 0.41 ), Pssbar( 0.18 ); // Broken SU(3) symmetry
681 G4cout<<"TarMinDiffMass "<< GetTar << 669 SetQuarkProbabilitiesAtGluonSplitUp( Puubar, Pddbar, Pssbar );
682 G4cout<<"TarMinNonDiffMass "<< GetTar <<
683 G4cout<<"AveragePt2 "<< GetAve <<
684 G4cout<<"ProbLogDistrPrD "<< GetPro <<
685 G4cout<<"ProbLogDistrTrD "<< GetPro <<
686 #endif <<
687 670
688 // Set parameters of nuclear destruction 671 // Set parameters of nuclear destruction
689 <<
690 if ( ProjectileabsPDGcode < 1000 ) { // Mes 672 if ( ProjectileabsPDGcode < 1000 ) { // Meson projectile
691 <<
692 const G4int indexTune = G4FTFTunings::Inst <<
693 <<
694 SetMaxNumberOfCollisions( Plab, 2.0 ); // 673 SetMaxNumberOfCollisions( Plab, 2.0 ); // 3.0 )
695 // << 674 //AR-18May2016 SetCofNuclearDestruction( 0.00481*G4double(NumberOfTargetNucleons)* // Uzhi 3.05.2015
696 // target destruction << 675 SetCofNuclearDestruction( 1.0* // AR-18May2016
697 // <<
698 /* original code ---> <<
699 SetCofNuclearDestruction( 0.00481*G4double <<
700 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + 676 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
701 <<
702 SetR2ofNuclearDestruction( 1.5*fermi*fermi 677 SetR2ofNuclearDestruction( 1.5*fermi*fermi );
703 SetDofNuclearDestruction( 0.3 ); 678 SetDofNuclearDestruction( 0.3 );
704 SetPt2ofNuclearDestruction( ( 0.035 + 0.04 679 SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
705 ( 1.0 680 ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
706 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV 681 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
707 SetExcitationEnergyPerWoundedNucleon( 40.0 682 SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
708 */ << 683 } else if ( ProjectilePDGcode < -1000 ) { // for anti-baryon projectile
709 double coeff = fArrayParCollMesonProj[inde << 684 SetMaxNumberOfCollisions( Plab, 2.0 ); // 3.0 )
710 // << 685 //AR-18May2016 SetCofNuclearDestruction( 0.00481*G4double(NumberOfTargetNucleons)* // Uzhi 3.05.2015
711 // NOTE (JVY): Set this switch to false/tr << 686 SetCofNuclearDestruction( 1.0* // AR-18May2016
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 687 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
742 SetR2ofNuclearDestruction( 1.5*fermi*fermi 688 SetR2ofNuclearDestruction( 1.5*fermi*fermi );
743 SetDofNuclearDestruction( 0.3 ); 689 SetDofNuclearDestruction( 0.3 );
744 SetPt2ofNuclearDestruction( ( 0.035 + 0.04 690 SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
745 ( 1.0 691 ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
746 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV 692 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
747 SetExcitationEnergyPerWoundedNucleon( 40.0 693 SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
748 if ( Plab < 2.0 ) { // 2 GeV/c 694 if ( Plab < 2.0 ) { // 2 GeV/c
749 // For slow anti-baryon we have to garan 695 // For slow anti-baryon we have to garanty putting on mass-shell
750 SetCofNuclearDestruction( 0.0 ); 696 SetCofNuclearDestruction( 0.0 );
751 SetR2ofNuclearDestruction( 1.5*fermi*fer << 697 SetR2ofNuclearDestruction( 1.5*fermi*fermi );
752 <<
753 SetDofNuclearDestruction( 0.01 ); 698 SetDofNuclearDestruction( 0.01 );
754 SetPt2ofNuclearDestruction( 0.035*GeV*Ge 699 SetPt2ofNuclearDestruction( 0.035*GeV*GeV );
755 SetMaxPt2ofNuclearDestruction( 0.04*GeV* 700 SetMaxPt2ofNuclearDestruction( 0.04*GeV*GeV );
>> 701 //SetExcitationEnergyPerWoundedNucleon( 0.0 ); // ?????
756 } 702 }
757 <<
758 } else { // Projectile baryon assumed 703 } else { // Projectile baryon assumed
759 << 704 SetMaxNumberOfCollisions( Plab, 2.0 ); // 3.0 )
760 // Below, in the call to the G4FTFTunings: << 705 //AR-18May2016 SetCofNuclearDestructionPr( 0.00481*G4double(AbsProjectileBaryonNumber)* // Uzhi 3.05.2015
761 // as projectile, instead of the real one, << 706 SetCofNuclearDestructionPr( 1.0* // AR-18May2016
762 // destruction, we assume for this categor << 707 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
763 // as for "baryon". << 708 //AR-18May2016 SetCofNuclearDestruction( 0.00481*G4double(NumberOfTargetNucleons)* // Uzhi 3.05.2015
764 const G4int indexTune = G4FTFTunings::Inst << 709 SetCofNuclearDestruction( 1.0* // AR-18May2016
765 << 710 G4Exp( 4.0*(Ylab - 2.1) )/( 1.0 + G4Exp( 4.0*(Ylab - 2.1) ) ) );
766 // NOTE (JVY) FIXME !!! Will decide later << 711 SetR2ofNuclearDestruction( 1.5*fermi*fermi );
767 // << 712 SetDofNuclearDestruction( 0.3 );
768 SetMaxNumberOfCollisions( Plab, 2.0 ); << 713 SetPt2ofNuclearDestruction( ( 0.035 + 0.04*G4Exp( 4.0*(Ylab - 2.5) )/
769 << 714 ( 1.0 + G4Exp( 4.0*(Ylab - 2.5) ) ) )*GeV*GeV );
770 // projectile destruction - does NOT reall << 715 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV );
771 // << 716 SetExcitationEnergyPerWoundedNucleon( 40.0*MeV );
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 } 717 }
817 718
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* 719 //SetCofNuclearDestruction( 0.47*G4Exp( 2.0*(Ylab - 2.5) )/( 1.0 + G4Exp( 2.0*(Ylab - 2.5) ) ) );
828 //SetPt2ofNuclearDestruction( ( 0.035 + 0.1* 720 //SetPt2ofNuclearDestruction( ( 0.035 + 0.1*G4Exp( 4.0*(Ylab - 3.0) )/( 1.0 + G4Exp( 4.0*(Ylab - 3.0) ) ) )*GeV*GeV );
829 721
830 //SetMagQuarkExchange( 120.0 ); // 210 722 //SetMagQuarkExchange( 120.0 ); // 210.0 PipP
831 //SetSlopeQuarkExchange( 2.0 ); 723 //SetSlopeQuarkExchange( 2.0 );
832 //SetDeltaProbAtQuarkExchange( 0.6 ); 724 //SetDeltaProbAtQuarkExchange( 0.6 );
833 //SetProjMinDiffMass( 0.7 ); // GeV 725 //SetProjMinDiffMass( 0.7 ); // GeV 1.1
834 //SetProjMinNonDiffMass( 0.7 ); // GeV 726 //SetProjMinNonDiffMass( 0.7 ); // GeV
835 //SetProbabilityOfProjDiff( 0.0); // 0.8 727 //SetProbabilityOfProjDiff( 0.0); // 0.85*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
836 //SetTarMinDiffMass( 1.1 ); // GeV 728 //SetTarMinDiffMass( 1.1 ); // GeV
837 //SetTarMinNonDiffMass( 1.1 ); // GeV 729 //SetTarMinNonDiffMass( 1.1 ); // GeV
838 //SetProbabilityOfTarDiff( 0.0 ); // 0.8 730 //SetProbabilityOfTarDiff( 0.0 ); // 0.85*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.5 ) ); // 40/32 X-dif/X-inel
839 731
840 //SetAveragePt2( 0.0 ); // GeV 732 //SetAveragePt2( 0.0 ); // GeV^2 0.3
841 //------------------------------------ 733 //------------------------------------
842 //SetProbabilityOfElasticScatt( 1.0, 1.0); 734 //SetProbabilityOfElasticScatt( 1.0, 1.0); //(Xtotal, Xelastic);
843 //SetProbabilityOfProjDiff( 1.0*0.62*G4Pow:: 735 //SetProbabilityOfProjDiff( 1.0*0.62*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.51 ) ); // 0->1
844 //SetProbabilityOfTarDiff( 4.0*0.62*G4Pow::G 736 //SetProbabilityOfTarDiff( 4.0*0.62*G4Pow::GetInstance()->powA( s/GeV/GeV, -0.51 ) ); // 2->4
845 //SetAveragePt2( 0.3 ); 737 //SetAveragePt2( 0.3 ); // (0.15)
846 //SetAvaragePt2ofElasticScattering( 0.0 ); 738 //SetAvaragePt2ofElasticScattering( 0.0 );
847 739
848 //SetMaxNumberOfCollisions( Plab, 6.0 ); //( 740 //SetMaxNumberOfCollisions( Plab, 6.0 ); //(4.0*(Plab + 0.01), Plab); // 6.0 );
849 //SetAveragePt2( 0.15 ); 741 //SetAveragePt2( 0.15 );
850 //SetCofNuclearDestruction(-1.);//( 0.75 ); 742 //SetCofNuclearDestruction(-1.);//( 0.75 ); // (0.25)
851 //SetExcitationEnergyPerWoundedNucleon(0.);/ 743 //SetExcitationEnergyPerWoundedNucleon(0.);//( 30.0*MeV ); // (75.0*MeV)
852 //SetDofNuclearDestruction(0.);//( 0.2 ); // 744 //SetDofNuclearDestruction(0.);//( 0.2 ); //0.4 // 0.3 0.5
853 745
854 /* << 746 //SetPt2ofNuclearDestruction(0.);//(2.*0.075*GeV*GeV); //( 0.3*GeV*GeV ); // (0.168*GeV*GeV)
855 SetAveragePt2(0.3); << 747 //SetMaxNumberOfCollisions( Plab, 78.0 ); // 3.0 )
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 748
866 //SetRadiusOfHNinteractions2( Xtotal/pi/10.0 << 749 //G4cout << "Cnd " << GetCofNuclearDestruction() << G4endl;
867 //SetRadiusOfHNinteractions2( (Xtotal - Xela << 750 //G4cout << "Dnd " << GetDofNuclearDestruction() << G4endl;
868 //SetProbabilityOfElasticScatt( 1.0, 0.0); << 751 //G4cout << "Pt2 " << GetPt2ofNuclearDestruction()/GeV/GeV << G4endl;
869 /* << 752 //G4int Uzhi; G4cin >> Uzhi;
870 G4cout << "Pt2 " << GetAveragePt2()<<" "<<Ge << 753 }
871 G4cout << "Cnd " << GetCofNuclearDestruction <<
872 G4cout << "Dnd " << GetDofNuclearDestruction <<
873 G4cout << "Pt2 " << GetPt2ofNuclearDestructi <<
874 */ <<
875 <<
876 } <<
877 754
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 755
902 //============================================ 756 //============================================================================
903 757
904 G4double G4FTFParameters::GetProcProb( const G 758 G4double G4FTFParameters::GetProcProb( const G4int ProcN, const G4double y ) {
905 G4double Prob( 0.0 ); 759 G4double Prob( 0.0 );
906 if ( y < ProcParams[ProcN][6] ) { 760 if ( y < ProcParams[ProcN][6] ) {
907 Prob = ProcParams[ProcN][5]; 761 Prob = ProcParams[ProcN][5];
908 if (Prob < 0.) Prob=0.; << 762 if(Prob < 0.) Prob=0.;
909 return Prob; 763 return Prob;
910 } 764 }
911 Prob = ProcParams[ProcN][0] * G4Exp( -ProcPa 765 Prob = ProcParams[ProcN][0] * G4Exp( -ProcParams[ProcN][1]*y ) +
912 ProcParams[ProcN][2] * G4Exp( -ProcPa 766 ProcParams[ProcN][2] * G4Exp( -ProcParams[ProcN][3]*y ) +
913 ProcParams[ProcN][4]; 767 ProcParams[ProcN][4];
914 if (Prob < 0.) Prob=0.; << 768 if(Prob < 0.) Prob=0.;
915 return Prob; 769 return Prob;
916 } 770 }
917 <<
918 //============================================ <<
919 <<
920 G4FTFParameters::~G4FTFParameters() { <<
921 if ( StringMass ) delete StringMass; <<
922 } <<
923 <<
924 //============================================ <<
925 <<
926 void G4FTFParameters::Reset() <<
927 { <<
928 FTFhNcmsEnergy = 0.0; <<
929 FTFXtotal = 0.0; <<
930 FTFXelastic = 0.0; <<
931 FTFXinelastic = 0.0; <<
932 FTFXannihilation = 0.0; <<
933 ProbabilityOfAnnihilation = 0.0; <<
934 ProbabilityOfElasticScatt = 0.0; <<
935 RadiusOfHNinteractions2 = 0.0; <<
936 FTFSlope = 0.0; <<
937 AvaragePt2ofElasticScattering = 0.0; <<
938 FTFGamma0 = 0.0; <<
939 DeltaProbAtQuarkExchange = 0.0; <<
940 ProbOfSameQuarkExchange = 0.0; <<
941 ProjMinDiffMass = 0.0; <<
942 ProjMinNonDiffMass = 0.0; <<
943 ProbLogDistrPrD = 0.0; <<
944 TarMinDiffMass = 0.0; <<
945 TarMinNonDiffMass = 0.0; <<
946 AveragePt2 = 0.0; <<
947 ProbLogDistr = 0.0; <<
948 Pt2kink = 0.0; <<
949 MaxNumberOfCollisions = 0.0; <<
950 ProbOfInelInteraction = 0.0; <<
951 CofNuclearDestructionPr = 0.0; <<
952 CofNuclearDestruction = 0.0; <<
953 R2ofNuclearDestruction = 0.0; <<
954 ExcitationEnergyPerWoundedNucleon = 0.0; <<
955 DofNuclearDestruction = 0.0; <<
956 Pt2ofNuclearDestruction = 0.0; <<
957 MaxPt2ofNuclearDestruction = 0.0; <<
958 <<
959 for ( G4int i = 0; i < 4; i++ ) { <<
960 for ( G4int j = 0; j < 7; j++ ) { <<
961 ProcParams[i][j] = 0.0; <<
962 } <<
963 } <<
964 <<
965 return; <<
966 } <<
967 <<
968 771