| 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 // 27 //
28 // Physics model class G4hhElastic 28 // Physics model class G4hhElastic
29 // 29 //
30 // 30 //
31 // G4 Model: qQ hadron hadron elastic scatteri 31 // G4 Model: qQ hadron hadron elastic scattering with 4-momentum balance
32 // 32 //
33 // 02.05.2014 V. Grichine 1-st version 33 // 02.05.2014 V. Grichine 1-st version
34 // 34 //
35 35
36 #include "G4hhElastic.hh" 36 #include "G4hhElastic.hh"
37 #include "G4ParticleTable.hh" 37 #include "G4ParticleTable.hh"
38 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh"
39 #include "G4IonTable.hh" 39 #include "G4IonTable.hh"
40 #include "G4NucleiProperties.hh" 40 #include "G4NucleiProperties.hh"
41 41
42 #include "Randomize.hh" 42 #include "Randomize.hh"
43 #include "G4Integrator.hh" 43 #include "G4Integrator.hh"
44 #include "globals.hh" 44 #include "globals.hh"
45 #include "G4PhysicalConstants.hh" 45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh"
47 47
48 #include "G4Proton.hh" 48 #include "G4Proton.hh"
49 #include "G4Neutron.hh" 49 #include "G4Neutron.hh"
50 #include "G4PionPlus.hh" 50 #include "G4PionPlus.hh"
51 #include "G4PionMinus.hh" 51 #include "G4PionMinus.hh"
52 52
53 #include "G4Element.hh" 53 #include "G4Element.hh"
54 #include "G4ElementTable.hh" 54 #include "G4ElementTable.hh"
55 #include "G4PhysicsTable.hh" 55 #include "G4PhysicsTable.hh"
56 #include "G4PhysicsLogVector.hh" 56 #include "G4PhysicsLogVector.hh"
57 #include "G4PhysicsFreeVector.hh" 57 #include "G4PhysicsFreeVector.hh"
58 58
59 #include "G4HadronNucleonXsc.hh" 59 #include "G4HadronNucleonXsc.hh"
60 60
61 #include "G4Pow.hh" 61 #include "G4Pow.hh"
62 62
63 #include "G4HadronicParameters.hh" 63 #include "G4HadronicParameters.hh"
64 64
65 using namespace std; 65 using namespace std;
66 66
67 67
68 ////////////////////////////////////////////// 68 /////////////////////////////////////////////////////////////////////////
69 // 69 //
70 // Tracking constructor. Target is proton 70 // Tracking constructor. Target is proton
71 71
72 72
73 G4hhElastic::G4hhElastic() 73 G4hhElastic::G4hhElastic()
74 : G4HadronElastic("HadrHadrElastic") 74 : G4HadronElastic("HadrHadrElastic")
75 { 75 {
76 SetMinEnergy( 1.*GeV ); 76 SetMinEnergy( 1.*GeV );
77 SetMaxEnergy( G4HadronicParameters::Instance 77 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
78 verboseLevel = 0; 78 verboseLevel = 0;
79 lowEnergyRecoilLimit = 100.*keV; 79 lowEnergyRecoilLimit = 100.*keV;
80 lowEnergyLimitQ = 0.0*GeV; 80 lowEnergyLimitQ = 0.0*GeV;
81 lowEnergyLimitHE = 0.0*GeV; 81 lowEnergyLimitHE = 0.0*GeV;
82 lowestEnergyLimit= 0.0*keV; 82 lowestEnergyLimit= 0.0*keV;
83 plabLowLimit = 20.0*MeV; 83 plabLowLimit = 20.0*MeV;
84 84
85 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 85 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
86 fInTkin=0; 86 fInTkin=0;
87 theProton = G4Proton::Proton(); 87 theProton = G4Proton::Proton();
88 theNeutron = G4Neutron::Neutron(); 88 theNeutron = G4Neutron::Neutron();
89 thePionPlus = G4PionPlus::PionPlus(); 89 thePionPlus = G4PionPlus::PionPlus();
90 thePionMinus= G4PionMinus::PionMinus(); 90 thePionMinus= G4PionMinus::PionMinus();
91 91
92 fTarget = G4Proton::Proton(); 92 fTarget = G4Proton::Proton();
93 fProjectile = 0; 93 fProjectile = 0;
94 fHadrNuclXsc = new G4HadronNucleonXsc(); 94 fHadrNuclXsc = new G4HadronNucleonXsc();
95 95
96 fEnergyBin = 200; 96 fEnergyBin = 200;
97 fBinT = 514; // 514; // 500; // 200; 97 fBinT = 514; // 514; // 500; // 200;
98 98
99 fEnergyVector = new G4PhysicsLogVector( the 99 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
100 100
101 fTableT = 0; 101 fTableT = 0;
102 fOldTkin = 0.; 102 fOldTkin = 0.;
103 SetParameters(); 103 SetParameters();
104 104
105 Initialise(); 105 Initialise();
106 } 106 }
107 107
108 108
109 ////////////////////////////////////////////// 109 /////////////////////////////////////////////////////////////////////////
110 // 110 //
111 // test constructor 111 // test constructor
112 112
113 113
114 G4hhElastic::G4hhElastic( G4ParticleDefinition 114 G4hhElastic::G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile, G4double plab)
115 : G4HadronElastic("HadrHadrElastic") 115 : G4HadronElastic("HadrHadrElastic")
116 { 116 {
117 SetMinEnergy( 1.*GeV ); 117 SetMinEnergy( 1.*GeV );
118 SetMaxEnergy( G4HadronicParameters::Instance 118 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
119 verboseLevel = 0; 119 verboseLevel = 0;
120 lowEnergyRecoilLimit = 100.*keV; 120 lowEnergyRecoilLimit = 100.*keV;
121 lowEnergyLimitQ = 0.0*GeV; 121 lowEnergyLimitQ = 0.0*GeV;
122 lowEnergyLimitHE = 0.0*GeV; 122 lowEnergyLimitHE = 0.0*GeV;
123 lowestEnergyLimit = 0.0*keV; 123 lowestEnergyLimit = 0.0*keV;
124 plabLowLimit = 20.0*MeV; 124 plabLowLimit = 20.0*MeV;
125 125
126 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 126 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
127 fInTkin=0; 127 fInTkin=0;
128 theProton = G4Proton::Proton(); 128 theProton = G4Proton::Proton();
129 theNeutron = G4Neutron::Neutron(); 129 theNeutron = G4Neutron::Neutron();
130 thePionPlus = G4PionPlus::PionPlus(); 130 thePionPlus = G4PionPlus::PionPlus();
131 thePionMinus= G4PionMinus::PionMinus(); 131 thePionMinus= G4PionMinus::PionMinus();
132 132
133 fTarget = target; 133 fTarget = target;
134 fProjectile = projectile; 134 fProjectile = projectile;
135 fMassTarg = fTarget->GetPDGMass(); 135 fMassTarg = fTarget->GetPDGMass();
136 fMassProj = fProjectile->GetPDGMass(); 136 fMassProj = fProjectile->GetPDGMass();
137 fMassSum2 = (fMassTarg+fMassProj)*(fMassTa 137 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
138 fMassDif2 = (fMassTarg-fMassProj)*(fMassTa 138 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
139 fHadrNuclXsc = new G4HadronNucleonXsc(); 139 fHadrNuclXsc = new G4HadronNucleonXsc();
140 140
141 fEnergyBin = 200; 141 fEnergyBin = 200;
142 fBinT = 514; // 200; 142 fBinT = 514; // 200;
143 143
144 fEnergyVector = new G4PhysicsLogVector( the 144 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
145 fTableT = 0; 145 fTableT = 0;
146 fOldTkin = 0.; 146 fOldTkin = 0.;
147 147
148 148
149 SetParameters(); 149 SetParameters();
150 SetParametersCMS( plab); 150 SetParametersCMS( plab);
151 } 151 }
152 152
153 153
154 ////////////////////////////////////////////// 154 /////////////////////////////////////////////////////////////////////////
155 // 155 //
156 // constructor used for low mass diffraction 156 // constructor used for low mass diffraction
157 157
158 158
159 G4hhElastic::G4hhElastic( G4ParticleDefinition 159 G4hhElastic::G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile)
160 : G4HadronElastic("HadrHadrElastic") 160 : G4HadronElastic("HadrHadrElastic")
161 { 161 {
162 SetMinEnergy( 1.*GeV ); 162 SetMinEnergy( 1.*GeV );
163 SetMaxEnergy( G4HadronicParameters::Instance 163 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
164 verboseLevel = 0; 164 verboseLevel = 0;
165 lowEnergyRecoilLimit = 100.*keV; 165 lowEnergyRecoilLimit = 100.*keV;
166 lowEnergyLimitQ = 0.0*GeV; 166 lowEnergyLimitQ = 0.0*GeV;
167 lowEnergyLimitHE = 0.0*GeV; 167 lowEnergyLimitHE = 0.0*GeV;
168 lowestEnergyLimit= 0.0*keV; 168 lowestEnergyLimit= 0.0*keV;
169 plabLowLimit = 20.0*MeV; 169 plabLowLimit = 20.0*MeV;
170 170
171 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0; 171 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
172 fInTkin=0; 172 fInTkin=0;
173 173
174 fTarget = target; // later vmg 174 fTarget = target; // later vmg
175 fProjectile = projectile; 175 fProjectile = projectile;
176 theProton = G4Proton::Proton(); 176 theProton = G4Proton::Proton();
177 theNeutron = G4Neutron::Neutron(); 177 theNeutron = G4Neutron::Neutron();
178 thePionPlus = G4PionPlus::PionPlus(); 178 thePionPlus = G4PionPlus::PionPlus();
179 thePionMinus= G4PionMinus::PionMinus(); 179 thePionMinus= G4PionMinus::PionMinus();
180 180
181 fTarget = G4Proton::Proton(); // later vmg 181 fTarget = G4Proton::Proton(); // later vmg
182 fMassTarg = fTarget->GetPDGMass(); 182 fMassTarg = fTarget->GetPDGMass();
183 fMassProj = fProjectile->GetPDGMass(); 183 fMassProj = fProjectile->GetPDGMass();
184 fMassSum2 = (fMassTarg+fMassProj)*(fMassTa 184 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
185 fMassDif2 = (fMassTarg-fMassProj)*(fMassTa 185 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
186 fHadrNuclXsc = new G4HadronNucleonXsc(); 186 fHadrNuclXsc = new G4HadronNucleonXsc();
187 187
188 fEnergyBin = 200; 188 fEnergyBin = 200;
189 fBinT = 514; // 514; // 500; // 200; 189 fBinT = 514; // 514; // 500; // 200;
190 190
191 fEnergyVector = new G4PhysicsLogVector( the 191 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
192 192
193 fTableT = 0; 193 fTableT = 0;
194 fOldTkin = 0.; 194 fOldTkin = 0.;
195 195
196 SetParameters(); 196 SetParameters();
197 } 197 }
198 198
199 199
200 200
201 ////////////////////////////////////////////// 201 //////////////////////////////////////////////////////////////////////////////
202 // 202 //
203 // Destructor 203 // Destructor
204 204
205 G4hhElastic::~G4hhElastic() 205 G4hhElastic::~G4hhElastic()
206 { 206 {
207 if ( fEnergyVector ) { 207 if ( fEnergyVector ) {
208 delete fEnergyVector; 208 delete fEnergyVector;
209 fEnergyVector = 0; 209 fEnergyVector = 0;
210 } 210 }
211 211
212 for ( std::vector<G4PhysicsTable*>::iterator 212 for ( std::vector<G4PhysicsTable*>::iterator it = fBankT.begin();
213 it != fBankT.end(); ++it ) { 213 it != fBankT.end(); ++it ) {
214 if ( (*it) ) (*it)->clearAndDestroy(); 214 if ( (*it) ) (*it)->clearAndDestroy();
215 delete *it; 215 delete *it;
216 *it = 0; 216 *it = 0;
217 } 217 }
218 fTableT = 0; 218 fTableT = 0;
219 if(fHadrNuclXsc) delete fHadrNuclXsc; 219 if(fHadrNuclXsc) delete fHadrNuclXsc;
220 } 220 }
221 221
222 ////////////////////////////////////////////// 222 /////////////////////////////////////////////////////////////////////////////
223 ///////////////////// Table preparation and r 223 ///////////////////// Table preparation and reading ////////////////////////
224 224
225 225
226 ////////////////////////////////////////////// 226 //////////////////////////////////////////////////////////////////////////////
227 // 227 //
228 // Initialisation for given particle on the pr 228 // Initialisation for given particle on the proton target
229 229
230 void G4hhElastic::Initialise() 230 void G4hhElastic::Initialise()
231 { 231 {
232 // pp,pn 232 // pp,pn
233 233
234 fProjectile = G4Proton::Proton(); 234 fProjectile = G4Proton::Proton();
235 BuildTableT(fTarget, fProjectile); 235 BuildTableT(fTarget, fProjectile);
236 fBankT.push_back(fTableT); // 0 236 fBankT.push_back(fTableT); // 0
237 237
238 // pi+-p 238 // pi+-p
239 239
240 fProjectile = G4PionPlus::PionPlus(); 240 fProjectile = G4PionPlus::PionPlus();
241 BuildTableT(fTarget, fProjectile); 241 BuildTableT(fTarget, fProjectile);
242 fBankT.push_back(fTableT); // 1 242 fBankT.push_back(fTableT); // 1
243 //K+-p 243 //K+-p
244 fProjectile = G4KaonPlus::KaonPlus(); 244 fProjectile = G4KaonPlus::KaonPlus();
245 BuildTableT(fTarget, fProjectile); 245 BuildTableT(fTarget, fProjectile);
246 fBankT.push_back(fTableT); // 2 246 fBankT.push_back(fTableT); // 2
247 247
248 } 248 }
249 249
250 ////////////////////////////////////////////// 250 ///////////////////////////////////////////////////////////////////////////////
251 // 251 //
252 // Build for given particle and proton table o 252 // Build for given particle and proton table of momentum transfers.
253 253
254 void G4hhElastic::BuildTableT( G4ParticleDefin 254 void G4hhElastic::BuildTableT( G4ParticleDefinition* target, G4ParticleDefinition* projectile) // , G4double plab)
255 { 255 {
256 G4int iTkin, jTransfer; 256 G4int iTkin, jTransfer;
257 G4double plab, Tkin, tMax; 257 G4double plab, Tkin, tMax;
258 G4double t1, t2, dt, delta = 0., sum = 0.; 258 G4double t1, t2, dt, delta = 0., sum = 0.;
259 259
260 fTarget = target; 260 fTarget = target;
261 fProjectile = projectile; 261 fProjectile = projectile;
262 fMassTarg = fTarget->GetPDGMass(); 262 fMassTarg = fTarget->GetPDGMass();
263 fMassProj = fProjectile->GetPDGMass(); 263 fMassProj = fProjectile->GetPDGMass();
264 fMassSum2 = (fMassTarg+fMassProj)*(fMassTa 264 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
265 fMassDif2 = (fMassTarg-fMassProj)*(fMassTa 265 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
266 266
267 G4Integrator<G4hhElastic,G4double(G4hhElasti 267 G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral;
268 // G4HadronNucleonXsc* hnXsc = new 268 // G4HadronNucleonXsc* hnXsc = new G4HadronNucleonXsc();
269 fTableT = new G4PhysicsTable(fEnergyBin); 269 fTableT = new G4PhysicsTable(fEnergyBin);
270 270
271 for( iTkin = 0; iTkin < fEnergyBin; iTkin++) 271 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
272 { 272 {
273 Tkin = fEnergyVector->GetLowEdgeEnergy(iT 273 Tkin = fEnergyVector->GetLowEdgeEnergy(iTkin);
274 plab = std::sqrt( Tkin*( Tkin + 2*fMassPr 274 plab = std::sqrt( Tkin*( Tkin + 2*fMassProj ) );
275 // G4DynamicParticle* theDynamicParticle 275 // G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(projectile,
276 // 276 // G4ParticleMomentum(0.,0.,1.),
277 // 277 // Tkin);
278 // fSigmaTot = fHadrNuclXsc->GetHadronNucl 278 // fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, target );
279 279
280 SetParametersCMS( plab ); 280 SetParametersCMS( plab );
281 281
282 tMax = 4.*fPcms*fPcms; 282 tMax = 4.*fPcms*fPcms;
283 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*Ge 283 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???
284 284
285 G4PhysicsFreeVector* vectorT = new G4Physi 285 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);
286 sum = 0.; 286 sum = 0.;
287 dt = tMax/fBinT; 287 dt = tMax/fBinT;
288 288
289 // for(j = 1; j < fBinT; j++) 289 // for(j = 1; j < fBinT; j++)
290 290
291 for( jTransfer = fBinT-1; jTransfer >= 1; 291 for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer--)
292 { 292 {
293 t1 = dt*(jTransfer-1); 293 t1 = dt*(jTransfer-1);
294 t2 = t1 + dt; 294 t2 = t1 + dt;
295 295
296 if( fMassProj > 900.*MeV ) // pp, pn 296 if( fMassProj > 900.*MeV ) // pp, pn
297 { 297 {
298 delta = integral.Legendre10(this, &G4h 298 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
299 // delta = integral.Legendre96(this, & 299 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
300 } 300 }
301 else // pi+-p, K+-p 301 else // pi+-p, K+-p
302 { 302 {
303 delta = integral.Legendre10(this, &G4h 303 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
304 // delta = integral.Legendre96(this, & 304 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
305 } 305 }
306 sum += delta; 306 sum += delta;
307 vectorT->PutValue( jTransfer-1, t1, sum 307 vectorT->PutValue( jTransfer-1, t1, sum ); // t2
308 } 308 }
309 // vectorT->PutValue( fBinT-1, dt*(fBinT-1 309 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
310 fTableT->insertAt( iTkin, vectorT ); 310 fTableT->insertAt( iTkin, vectorT );
311 // delete theDynamicParticle; 311 // delete theDynamicParticle;
312 } 312 }
313 // delete hnXsc; 313 // delete hnXsc;
314 314
315 return; 315 return;
316 } 316 }
317 317
318 ////////////////////////////////////////////// 318 ////////////////////////////////////////////////////////////////////////////
319 // 319 //
320 // Return inv momentum transfer -t > 0 from in 320 // Return inv momentum transfer -t > 0 from initialisation table
321 321
322 G4double G4hhElastic::SampleInvariantT( const 322 G4double G4hhElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p,
323 323 G4int, G4int )
324 { 324 {
325 G4int iTkin, iTransfer; 325 G4int iTkin, iTransfer;
326 G4double t, t2, position, m1 = aParticle->Ge 326 G4double t, t2, position, m1 = aParticle->GetPDGMass();
327 G4double Tkin = std::sqrt(m1*m1+p*p) - m1; 327 G4double Tkin = std::sqrt(m1*m1+p*p) - m1;
328 328
329 if( aParticle == G4Proton::Proton() || aPart 329 if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
330 { 330 {
331 fTableT = fBankT[0]; 331 fTableT = fBankT[0];
332 } 332 }
333 if( aParticle == G4PionPlus::PionPlus() || a 333 if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
334 { 334 {
335 fTableT = fBankT[1]; 335 fTableT = fBankT[1];
336 } 336 }
337 if( aParticle == G4KaonPlus::KaonPlus() || a 337 if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
338 { 338 {
339 fTableT = fBankT[2]; 339 fTableT = fBankT[2];
340 } 340 }
341 341
342 G4double delta = std::abs(Tkin - fOldTkin 342 G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
343 G4double deltaMax = 1.e-2; 343 G4double deltaMax = 1.e-2;
344 344
345 if ( delta < deltaMax ) iTkin = fInTkin; 345 if ( delta < deltaMax ) iTkin = fInTkin;
346 else 346 else
347 { 347 {
348 for( iTkin = 0; iTkin < fEnergyBin; iTkin+ 348 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
349 { 349 {
350 if( Tkin < fEnergyVector->GetLowEdgeEner 350 if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
351 } 351 }
352 } 352 }
353 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi 353 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
354 if ( iTkin < 0 ) iTkin = 0; // aga 354 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
355 355
356 fOldTkin = Tkin; 356 fOldTkin = Tkin;
357 fInTkin = iTkin; 357 fInTkin = iTkin;
358 358
359 if (iTkin == fEnergyBin -1 || iTkin == 0 ) 359 if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges
360 { 360 {
361 position = (*(*fTableT)(iTkin))(0)*G4Unifo 361 position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
362 362
363 // G4cout<<"position = "<<position<<G4endl 363 // G4cout<<"position = "<<position<<G4endl;
364 364
365 for(iTransfer = 0; iTransfer < fBinT-1; iT 365 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
366 { 366 {
367 if( position >= (*(*fTableT)(iTkin))(iTr 367 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
368 } 368 }
369 if (iTransfer >= fBinT-1) iTransfer = fBin 369 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
370 370
371 // G4cout<<"iTransfer = "<<iTransfer<<G4en 371 // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
372 372
373 t = GetTransfer(iTkin, iTransfer, position 373 t = GetTransfer(iTkin, iTransfer, position);
374 374
375 // G4cout<<"t = "<<t<<G4endl; 375 // G4cout<<"t = "<<t<<G4endl;
376 } 376 }
377 else // Tkin inside between energy table ed 377 else // Tkin inside between energy table edges
378 { 378 {
379 // position = (*(*fTableT)(iTkin))(fBinT-2 379 // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();
380 position = (*(*fTableT)(iTkin))(0)*G4Unifo 380 position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
381 381
382 // G4cout<<"position = "<<position<<G4endl 382 // G4cout<<"position = "<<position<<G4endl;
383 383
384 for(iTransfer = 0; iTransfer < fBinT-1; iT 384 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
385 { 385 {
386 // if( position < (*(*fTableT)(iTkin))(i 386 // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
387 if( position >= (*(*fTableT)(iTkin))(iTr 387 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
388 } 388 }
389 if (iTransfer >= fBinT-1) iTransfer = fBin 389 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
390 390
391 // G4cout<<"iTransfer = "<<iTransfer<<G4en 391 // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
392 392
393 t2 = GetTransfer(iTkin, iTransfer, positi 393 t2 = GetTransfer(iTkin, iTransfer, position);
394 return t2; 394 return t2;
395 /* 395 /*
396 G4double t1, E1, E2, W, W1, W2; 396 G4double t1, E1, E2, W, W1, W2;
397 // G4cout<<"t2 = "<<t2<<G4endl; 397 // G4cout<<"t2 = "<<t2<<G4endl;
398 398
399 E2 = fEnergyVector->GetLowEdgeEnergy(iTkin 399 E2 = fEnergyVector->GetLowEdgeEnergy(iTkin);
400 400
401 // G4cout<<"E2 = "<<E2<<G4endl; 401 // G4cout<<"E2 = "<<E2<<G4endl;
402 402
403 iTkin--; 403 iTkin--;
404 404
405 // position = (*(*fTableT)(iTkin))(fBinT-2 405 // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();
406 406
407 // G4cout<<"position = "<<position<<G4endl 407 // G4cout<<"position = "<<position<<G4endl;
408 408
409 for(iTransfer = 0; iTransfer < fBinT-1; iT 409 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
410 { 410 {
411 // if( position < (*(*fTableT)(iTkin))(i 411 // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
412 if( position >= (*(*fTableT)(iTkin))(iTr 412 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
413 } 413 }
414 if (iTransfer >= fBinT-1) iTransfer = fBin 414 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
415 415
416 t1 = GetTransfer(iTkin, iTransfer, positi 416 t1 = GetTransfer(iTkin, iTransfer, position);
417 417
418 // G4cout<<"t1 = "<<t1<<G4endl; 418 // G4cout<<"t1 = "<<t1<<G4endl;
419 419
420 E1 = fEnergyVector->GetLowEdgeEnergy(iTkin 420 E1 = fEnergyVector->GetLowEdgeEnergy(iTkin);
421 421
422 // G4cout<<"E1 = "<<E1<<G4endl; 422 // G4cout<<"E1 = "<<E1<<G4endl;
423 423
424 W = 1.0/(E2 - E1); 424 W = 1.0/(E2 - E1);
425 W1 = (E2 - Tkin)*W; 425 W1 = (E2 - Tkin)*W;
426 W2 = (Tkin - E1)*W; 426 W2 = (Tkin - E1)*W;
427 427
428 t = W1*t1 + W2*t2; 428 t = W1*t1 + W2*t2;
429 */ 429 */
430 } 430 }
431 return t; 431 return t;
432 } 432 }
433 433
434 434
435 ////////////////////////////////////////////// 435 ////////////////////////////////////////////////////////////////////////////
436 // 436 //
437 // Return inv momentum transfer -t > 0 from in 437 // Return inv momentum transfer -t > 0 from initialisation table
438 438
439 G4double G4hhElastic::SampleBisectionalT( cons 439 G4double G4hhElastic::SampleBisectionalT( const G4ParticleDefinition* aParticle, G4double p)
440 { 440 {
441 G4int iTkin, iTransfer; 441 G4int iTkin, iTransfer;
442 G4double t, position, m1 = aParticle->GetPD 442 G4double t, position, m1 = aParticle->GetPDGMass();
443 G4double Tkin = std::sqrt(m1*m1+p*p) - m1; 443 G4double Tkin = std::sqrt(m1*m1+p*p) - m1;
444 444
445 if( aParticle == G4Proton::Proton() || aPart 445 if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
446 { 446 {
447 fTableT = fBankT[0]; 447 fTableT = fBankT[0];
448 } 448 }
449 if( aParticle == G4PionPlus::PionPlus() || a 449 if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
450 { 450 {
451 fTableT = fBankT[1]; 451 fTableT = fBankT[1];
452 } 452 }
453 if( aParticle == G4KaonPlus::KaonPlus() || a 453 if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
454 { 454 {
455 fTableT = fBankT[2]; 455 fTableT = fBankT[2];
456 } 456 }
457 G4double delta = std::abs(Tkin - fOldTkin 457 G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
458 G4double deltaMax = 1.e-2; 458 G4double deltaMax = 1.e-2;
459 459
460 if ( delta < deltaMax ) iTkin = fInTkin; 460 if ( delta < deltaMax ) iTkin = fInTkin;
461 else 461 else
462 { 462 {
463 for( iTkin = 0; iTkin < fEnergyBin; iTkin+ 463 for( iTkin = 0; iTkin < fEnergyBin; iTkin++ )
464 { 464 {
465 if( Tkin < fEnergyVector->GetLowEdgeEner 465 if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
466 } 466 }
467 } 467 }
468 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi 468 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
469 if ( iTkin < 0 ) iTkin = 0; // aga 469 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
470 470
471 fOldTkin = Tkin; 471 fOldTkin = Tkin;
472 fInTkin = iTkin; 472 fInTkin = iTkin;
473 473
474 if (iTkin == fEnergyBin -1 || iTkin == 0 ) 474 if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges
475 { 475 {
476 position = (*(*fTableT)(iTkin))(0)*G4Unifo 476 position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
477 477
478 for(iTransfer = 0; iTransfer < fBinT-1; iT 478 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
479 { 479 {
480 if( position >= (*(*fTableT)(iTkin))(iTr 480 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
481 } 481 }
482 if (iTransfer >= fBinT-1) iTransfer = fBin 482 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
483 483
484 t = GetTransfer(iTkin, iTransfer, position 484 t = GetTransfer(iTkin, iTransfer, position);
485 485
486 486
487 } 487 }
488 else // Tkin inside between energy table ed 488 else // Tkin inside between energy table edges
489 { 489 {
490 G4double rand = G4UniformRand(); 490 G4double rand = G4UniformRand();
491 position = (*(*fTableT)(iTkin))(0)*rand; 491 position = (*(*fTableT)(iTkin))(0)*rand;
492 492
493 // 493 //
494 // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBi 494 // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBinT-2);
495 G4int sTransfer = 0, fTransfer = fBinT - 495 G4int sTransfer = 0, fTransfer = fBinT - 2, dTransfer = fTransfer - sTransfer;
496 G4double y2; 496 G4double y2;
497 497
498 for( iTransfer = 0; iTransfer < fBinT - 1; 498 for( iTransfer = 0; iTransfer < fBinT - 1; iTransfer++ )
499 { 499 {
500 // dTransfer %= 2; 500 // dTransfer %= 2;
501 dTransfer /= 2; 501 dTransfer /= 2;
502 // dTransfer *= 0.5; 502 // dTransfer *= 0.5;
503 y2 = (*(*fTableT)(iTkin))( sTransfer + d 503 y2 = (*(*fTableT)(iTkin))( sTransfer + dTransfer );
504 504
505 if( y2 > position ) sTransfer += dTransf 505 if( y2 > position ) sTransfer += dTransfer;
506 506
507 // if( dTransfer <= 1 ) break; 507 // if( dTransfer <= 1 ) break;
508 if( dTransfer < 1 ) break; 508 if( dTransfer < 1 ) break;
509 } 509 }
510 t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sT 510 t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sTransfer); // +(-0.5+rand)*(*fTableT)(iTkin)->GetLowEdgeEnergy(3);
511 } 511 }
512 return t; 512 return t;
513 } 513 }
514 514
515 515
516 ////////////////////////////////////////////// 516 ///////////////////////////////////////////////////////////////////////////////
517 // 517 //
518 // Build for given particle and proton table o 518 // Build for given particle and proton table of momentum transfers.
519 519
520 void G4hhElastic::BuildTableTest( G4ParticleDe 520 void G4hhElastic::BuildTableTest( G4ParticleDefinition* target, G4ParticleDefinition* projectile, G4double plab)
521 { 521 {
522 G4int jTransfer; 522 G4int jTransfer;
523 G4double tMax; // , sQq, sQG; 523 G4double tMax; // , sQq, sQG;
524 G4double t1, t2, dt, delta = 0., sum = 0. ; 524 G4double t1, t2, dt, delta = 0., sum = 0. ; // , threshold;
525 525
526 fTarget = target; 526 fTarget = target;
527 fProjectile = projectile; 527 fProjectile = projectile;
528 fMassTarg = fTarget->GetPDGMass(); 528 fMassTarg = fTarget->GetPDGMass();
529 fMassProj = fProjectile->GetPDGMass(); 529 fMassProj = fProjectile->GetPDGMass();
530 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg 530 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
531 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg 531 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
532 fSpp = fMassProj*fMassProj + fMassTarg*fMas 532 fSpp = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
533 fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp 533 fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
534 534
535 G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fM 535 G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fMassProj = "<<fMassProj<<" MeV"<<G4endl;
536 tMax = 4.*fPcms*fPcms; 536 tMax = 4.*fPcms*fPcms;
537 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; 537 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???
538 538
539 539
540 G4Integrator<G4hhElastic,G4double(G4hhElasti 540 G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral;
541 fTableT = new G4PhysicsTable(1); 541 fTableT = new G4PhysicsTable(1);
542 G4PhysicsFreeVector* vectorT = new G4Physics 542 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);
543 543
544 sum = 0.; 544 sum = 0.;
545 dt = tMax/G4double(fBinT); 545 dt = tMax/G4double(fBinT);
546 G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV; 546 G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV; fPcms = "<<fPcms/GeV
547 <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt 547 <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt = "<<dt/GeV/GeV<<" GeV2"<<G4endl;
548 548
549 // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fR 549 // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fRB*GeV<<G4endl;
550 550
551 // for(jTransfer = 1; jTransfer < fBinT; j 551 // for(jTransfer = 1; jTransfer < fBinT; jTransfer++)
552 for( jTransfer = fBinT-1; jTransfer >= 1; jT 552 for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer-- )
553 { 553 {
554 t1 = dt*(jTransfer-1); 554 t1 = dt*(jTransfer-1);
555 t2 = t1 + dt; 555 t2 = t1 + dt;
556 556
557 if( fMassProj > 900.*MeV ) // pp, pn 557 if( fMassProj > 900.*MeV ) // pp, pn
558 { 558 {
559 delta = integral.Legendre10(this, &G4hhE 559 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
560 // threshold = integral.Legendre96(this, 560 // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, tMax);
561 } 561 }
562 else // pi+-p, K+-p 562 else // pi+-p, K+-p
563 { 563 {
564 delta = integral.Legendre10(this, &G4hhE 564 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
565 // threshold = integral.Legendre96(this, 565 // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, tMax);
566 // delta = integral.Legendre96(this, &G4 566 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
567 } 567 }
568 sum += delta; 568 sum += delta;
569 // G4cout<<delta<<"\t"<<sum<<"\t"<<thresho 569 // G4cout<<delta<<"\t"<<sum<<"\t"<<threshold<<G4endl;
570 570
571 // sQq = GetdsdtF123(q1); 571 // sQq = GetdsdtF123(q1);
572 // sQG = GetdsdtF123qQgG(q1); 572 // sQG = GetdsdtF123qQgG(q1);
573 // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/milli 573 // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/millibarn<<"\t"<<sQq*GeV*GeV/millibarn<<G4endl;
574 // G4cout<<"sum = "<<sum<<", "; 574 // G4cout<<"sum = "<<sum<<", ";
575 575
576 vectorT->PutValue( jTransfer-1, t1, sum ); 576 vectorT->PutValue( jTransfer-1, t1, sum ); // t2
577 } 577 }
578 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 578 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
579 fTableT->insertAt( 0, vectorT ); 579 fTableT->insertAt( 0, vectorT );
580 fBankT.push_back( fTableT ); // 0 580 fBankT.push_back( fTableT ); // 0
581 581
582 // for(jTransfer = 0; jTransfer < fBinT-1; j 582 // for(jTransfer = 0; jTransfer < fBinT-1; jTransfer++)
583 // G4cout<<(*(*fTableT)(0))(jTransfer)/sum 583 // G4cout<<(*(*fTableT)(0))(jTransfer)/sum<<"\t\t"<<G4Pow::GetInstance()->powN(2.,-jTransfer)<<G4endl;
584 584
585 return; 585 return;
586 } 586 }
587 587
588 588
589 ////////////////////////////////////////////// 589 ////////////////////////////////////////////////////////////////////////////
590 // 590 //
591 // Return inv momentum transfer -t > 0 from in 591 // Return inv momentum transfer -t > 0 from initialisation table
592 592
593 G4double G4hhElastic::SampleTest(G4double tMin 593 G4double G4hhElastic::SampleTest(G4double tMin ) // const G4ParticleDefinition* aParticle, )
594 { 594 {
595 G4int iTkin, iTransfer, iTmin; 595 G4int iTkin, iTransfer, iTmin;
596 G4double t, position; 596 G4double t, position;
597 // G4double qMin = std::sqrt(tMin); 597 // G4double qMin = std::sqrt(tMin);
598 598
599 fTableT = fBankT[0]; 599 fTableT = fBankT[0];
600 iTkin = 0; 600 iTkin = 0;
601 601
602 for(iTransfer = 0; iTransfer < fBinT-1; iTra 602 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
603 { 603 {
604 // if( qMin <= (*fTableT)(iTkin)->GetLowEd 604 // if( qMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
605 if( tMin <= (*fTableT)(iTkin)->GetLowEdgeE 605 if( tMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
606 } 606 }
607 iTmin = iTransfer-1; 607 iTmin = iTransfer-1;
608 if(iTmin < 0 ) iTmin = 0; 608 if(iTmin < 0 ) iTmin = 0;
609 609
610 position = (*(*fTableT)(iTkin))(iTmin)*G4Uni 610 position = (*(*fTableT)(iTkin))(iTmin)*G4UniformRand();
611 611
612 for( iTmin = 0; iTransfer < fBinT-1; iTransf 612 for( iTmin = 0; iTransfer < fBinT-1; iTransfer++)
613 { 613 {
614 if( position > (*(*fTableT)(iTkin))(iTrans 614 if( position > (*(*fTableT)(iTkin))(iTransfer) ) break;
615 } 615 }
616 if (iTransfer >= fBinT-1) iTransfer = fBinT- 616 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
617 617
618 t = GetTransfer(iTkin, iTransfer, position); 618 t = GetTransfer(iTkin, iTransfer, position);
619 619
620 return t; 620 return t;
621 } 621 }
622 622
623 623
624 ////////////////////////////////////////////// 624 /////////////////////////////////////////////////////////////////////////////////
625 // 625 //
626 // Check with PAI sampling 626 // Check with PAI sampling
627 627
628 G4double 628 G4double
629 G4hhElastic:: GetTransfer( G4int iTkin, G4int 629 G4hhElastic:: GetTransfer( G4int iTkin, G4int iTransfer, G4double position )
630 { 630 {
631 G4double x1, x2, y1, y2, randTransfer, delta 631 G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;
632 632
633 if( iTransfer == 0 ) 633 if( iTransfer == 0 )
634 { 634 {
635 randTransfer = (*fTableT)(iTkin)->GetLowEd 635 randTransfer = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);
636 // iTransfer++; 636 // iTransfer++;
637 } 637 }
638 else 638 else
639 { 639 {
640 if ( iTransfer >= G4int((*fTableT)(iTkin)- 640 if ( iTransfer >= G4int((*fTableT)(iTkin)->GetVectorLength()) )
641 { 641 {
642 iTransfer = G4int((*fTableT)(iTkin)->Get << 642 iTransfer = (*fTableT)(iTkin)->GetVectorLength() - 1;
643 } 643 }
644 y1 = (*(*fTableT)(iTkin))(iTransfer-1); 644 y1 = (*(*fTableT)(iTkin))(iTransfer-1);
645 y2 = (*(*fTableT)(iTkin))(iTransfer); 645 y2 = (*(*fTableT)(iTkin))(iTransfer);
646 646
647 x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(i 647 x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer-1);
648 x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(i 648 x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);
649 649
650 delta = y2 - y1; 650 delta = y2 - y1;
651 mean = y2 + y1; 651 mean = y2 + y1;
652 652
653 if ( x1 == x2 ) randTransfer = x2; 653 if ( x1 == x2 ) randTransfer = x2;
654 else 654 else
655 { 655 {
656 // if ( y1 == y2 ) 656 // if ( y1 == y2 )
657 if ( delta < epsilon*mean ) 657 if ( delta < epsilon*mean )
658 randTransfer = x1 + ( x2 - x1 )*G4U 658 randTransfer = x1 + ( x2 - x1 )*G4UniformRand();
659 else randTransfer = x1 + ( position - y1 659 else randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
660 } 660 }
661 } 661 }
662 return randTransfer; 662 return randTransfer;
663 } 663 }
664 664
665 const G4double G4hhElastic::theNuclNuclData[19 665 const G4double G4hhElastic::theNuclNuclData[19][6] =
666 { 666 {
667 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1 667 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1/GeV, fBq, fBQ, fImCof
668 668
669 { 2.76754, 4.8, 4.8, 0.05, 0.742441, 10.5 669 { 2.76754, 4.8, 4.8, 0.05, 0.742441, 10.5 }, // pp 3GeV/c
670 { 3.07744, 5.4, 5.4, 0.02, 0.83818, 6.5 670 { 3.07744, 5.4, 5.4, 0.02, 0.83818, 6.5 }, // pp 4GeV/c
671 { 3.36305, 5.2, 5.2, 0.02, 0.838893, 7.5 671 { 3.36305, 5.2, 5.2, 0.02, 0.838893, 7.5 }, // np 5GeV/c
672 { 4.32941, 6, 6, 0.03, 0.769389, 7.5 }, 672 { 4.32941, 6, 6, 0.03, 0.769389, 7.5 }, // np 9 GeV/c
673 { 4.62126, 6, 6, 0.03, 0.770111, 6.5 }, 673 { 4.62126, 6, 6, 0.03, 0.770111, 6.5 }, // pp 10.4 GeV/c
674 674
675 { 5.47416, 4.5, 4.5, 0.03, 0.813185, 7.5 675 { 5.47416, 4.5, 4.5, 0.03, 0.813185, 7.5 }, // np 15 GeV/c
676 { 6.15088, 6.5, 6.5, 0.02, 0.799539, 6.5 676 { 6.15088, 6.5, 6.5, 0.02, 0.799539, 6.5 }, // pp 19.2 GeV/c
677 { 6.77474, 5.2, 5.2, 0.03, 0.784901, 7.5 677 { 6.77474, 5.2, 5.2, 0.03, 0.784901, 7.5 }, // np 23.5 GeV/c
678 { 9.77775, 7, 7, 0.03, 0.742531, 6.5 }, 678 { 9.77775, 7, 7, 0.03, 0.742531, 6.5 }, // pp 50 GeV/c
679 // {9.77775, 7, 7, 0.011, 679 // {9.77775, 7, 7, 0.011, 0.84419, 4.5 }, // pp 50 GeV/c
680 { 10.4728, 5.2, 5.2, 0.03, 0.780439, 7.5 680 { 10.4728, 5.2, 5.2, 0.03, 0.780439, 7.5 }, // np 57.5 GeV/c
681 681
682 { 13.7631, 7, 7, 0.008, 0.8664, 682 { 13.7631, 7, 7, 0.008, 0.8664, 5.0 }, // pp 100 GeV/c
683 { 19.4184, 6.8, 6.8, 0.009, 0.861337, 2. 683 { 19.4184, 6.8, 6.8, 0.009, 0.861337, 2.5 }, // pp 200 GeV/c
684 { 23.5, 6.8, 6.8, 0.007, 0.878112, 1.5 } 684 { 23.5, 6.8, 6.8, 0.007, 0.878112, 1.5 }, // pp 23.5 GeV
685 // {24.1362, 6.4, 6.4, 0.09 685 // {24.1362, 6.4, 6.4, 0.09, 0.576215, 7.5 }, // np 309.5 GeV/c
686 { 24.1362, 7.2, 7.2, 0.008, 0.864745, 5. 686 { 24.1362, 7.2, 7.2, 0.008, 0.864745, 5.5 },
687 { 52.8, 6.8, 6.8, 0.008, 0.871929, 1.5 } 687 { 52.8, 6.8, 6.8, 0.008, 0.871929, 1.5 }, // pp 58.2 GeV
688 688
689 { 546, 7.4, 7.4, 0.013, 0.845877, 689 { 546, 7.4, 7.4, 0.013, 0.845877, 5.5 }, // pb-p 546 GeV
690 { 1960, 7.8, 7.8, 0.022, 0.809062, 7.5 } 690 { 1960, 7.8, 7.8, 0.022, 0.809062, 7.5 }, // pb-p 1960 GeV
691 { 7000, 8, 8, 0.024, 0.820441, 5.5 }, / 691 { 7000, 8, 8, 0.024, 0.820441, 5.5 }, // pp TOTEM 7 TeV
692 { 13000, 8.5, 8.5, 0.03, 0.796721, 10.5 692 { 13000, 8.5, 8.5, 0.03, 0.796721, 10.5 } // pp TOTEM 13 TeV
693 693
694 }; 694 };
695 695
696 ////////////////////////////////////////////// 696 //////////////////////////////////////////////////////////////////////////////////
697 697
698 const G4double G4hhElastic::thePiKaNuclData[8] 698 const G4double G4hhElastic::thePiKaNuclData[8][6] =
699 { 699 {
700 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1 700 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1/GeV, fBq, fBQ, fImCof
701 701
702 { 2.5627, 3.8, 3.3, 0.22, 0.222, 702 { 2.5627, 3.8, 3.3, 0.22, 0.222, 1.5 }, // pipp 3.017 GeV/c
703 { 2.93928, 4.3, 3.8, 0.2, 0.2506 703 { 2.93928, 4.3, 3.8, 0.2, 0.250601, 1.3 }, // pipp 4.122 GeV/c
704 { 3.22326, 4.8, 4.3, 0.13, 0.32751, 2.5 704 { 3.22326, 4.8, 4.3, 0.13, 0.32751, 2.5 }, // pipp 5.055 GeV/c
705 { 7.80704, 5.5, 5, 0.13, 0.340631, 2.5 }, 705 { 7.80704, 5.5, 5, 0.13, 0.340631, 2.5 }, // pipp 32 GeV/c
706 { 9.7328, 5, 4.5, 0.05, 0.416319, 5.5 }, 706 { 9.7328, 5, 4.5, 0.05, 0.416319, 5.5 }, // pipp 50 GeV/c
707 707
708 { 13.7315, 5.3, 4.8, 0.05, 0.418426, 5.5 708 { 13.7315, 5.3, 4.8, 0.05, 0.418426, 5.5 }, // pipp 100 GeV/c
709 { 16.6359, 6.3, 5.8, 0.05, 0.423817, 5.5 709 { 16.6359, 6.3, 5.8, 0.05, 0.423817, 5.5 }, // pipp 147 GeV/c
710 { 19.3961, 5, 4.5, 0.05, 0.413477, 3.5 } 710 { 19.3961, 5, 4.5, 0.05, 0.413477, 3.5 } // pimp 200 GeV/c
711 711
712 }; 712 };
713 713
714 // 714 //
715 // 715 //
716 ////////////////////////////////////////////// 716 /////////////////////////////////////////////////////////////////////////////////
717 717