Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4hhElastic.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27//
28// Physics model class G4hhElastic
29//
30//
31// G4 Model: qQ hadron hadron elastic scattering with 4-momentum balance
32//
33// 02.05.2014 V. Grichine 1-st version
34//
35
36#include "G4hhElastic.hh"
37#include "G4ParticleTable.hh"
38#include "G4ParticleDefinition.hh"
39#include "G4IonTable.hh"
40#include "G4NucleiProperties.hh"
41
42#include "Randomize.hh"
43#include "G4Integrator.hh"
44#include "globals.hh"
45#include "G4PhysicalConstants.hh"
46#include "G4SystemOfUnits.hh"
47
48#include "G4Proton.hh"
49#include "G4Neutron.hh"
50#include "G4PionPlus.hh"
51#include "G4PionMinus.hh"
52
53#include "G4Element.hh"
54#include "G4ElementTable.hh"
55#include "G4PhysicsTable.hh"
56#include "G4PhysicsLogVector.hh"
57#include "G4PhysicsFreeVector.hh"
58
59#include "G4HadronNucleonXsc.hh"
60
61#include "G4Pow.hh"
62
63#include "G4HadronicParameters.hh"
64
65using namespace std;
66
67
68/////////////////////////////////////////////////////////////////////////
69//
70// Tracking constructor. Target is proton
71
72
73G4hhElastic::G4hhElastic()
74 : G4HadronElastic("HadrHadrElastic")
75{
76 SetMinEnergy( 1.*GeV );
77 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
78 verboseLevel = 0;
79 lowEnergyRecoilLimit = 100.*keV;
80 lowEnergyLimitQ = 0.0*GeV;
81 lowEnergyLimitHE = 0.0*GeV;
82 lowestEnergyLimit= 0.0*keV;
83 plabLowLimit = 20.0*MeV;
84
85 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
86 fInTkin=0;
87 theProton = G4Proton::Proton();
88 theNeutron = G4Neutron::Neutron();
89 thePionPlus = G4PionPlus::PionPlus();
90 thePionMinus= G4PionMinus::PionMinus();
91
92 fTarget = G4Proton::Proton();
93 fProjectile = 0;
94 fHadrNuclXsc = new G4HadronNucleonXsc();
95
96 fEnergyBin = 200;
97 fBinT = 514; // 514; // 500; // 200;
98
99 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
100
101 fTableT = 0;
102 fOldTkin = 0.;
103 SetParameters();
104
105 Initialise();
106}
107
108
109/////////////////////////////////////////////////////////////////////////
110//
111// test constructor
112
113
114G4hhElastic::G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile, G4double plab)
115 : G4HadronElastic("HadrHadrElastic")
116{
117 SetMinEnergy( 1.*GeV );
118 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
119 verboseLevel = 0;
120 lowEnergyRecoilLimit = 100.*keV;
121 lowEnergyLimitQ = 0.0*GeV;
122 lowEnergyLimitHE = 0.0*GeV;
123 lowestEnergyLimit = 0.0*keV;
124 plabLowLimit = 20.0*MeV;
125
126 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
127 fInTkin=0;
128 theProton = G4Proton::Proton();
129 theNeutron = G4Neutron::Neutron();
130 thePionPlus = G4PionPlus::PionPlus();
131 thePionMinus= G4PionMinus::PionMinus();
132
133 fTarget = target;
134 fProjectile = projectile;
135 fMassTarg = fTarget->GetPDGMass();
136 fMassProj = fProjectile->GetPDGMass();
137 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
138 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
139 fHadrNuclXsc = new G4HadronNucleonXsc();
140
141 fEnergyBin = 200;
142 fBinT = 514; // 200;
143
144 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
145 fTableT = 0;
146 fOldTkin = 0.;
147
148
149 SetParameters();
150 SetParametersCMS( plab);
151}
152
153
154/////////////////////////////////////////////////////////////////////////
155//
156// constructor used for low mass diffraction
157
158
159G4hhElastic::G4hhElastic( G4ParticleDefinition* target, G4ParticleDefinition* projectile)
160 : G4HadronElastic("HadrHadrElastic")
161{
162 SetMinEnergy( 1.*GeV );
163 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
164 verboseLevel = 0;
165 lowEnergyRecoilLimit = 100.*keV;
166 lowEnergyLimitQ = 0.0*GeV;
167 lowEnergyLimitHE = 0.0*GeV;
168 lowestEnergyLimit= 0.0*keV;
169 plabLowLimit = 20.0*MeV;
170
171 fRhoReIm=fSigmaTot=fOptRatio=fSpp=fPcms=0.0;
172 fInTkin=0;
173
174 fTarget = target; // later vmg
175 fProjectile = projectile;
176 theProton = G4Proton::Proton();
177 theNeutron = G4Neutron::Neutron();
178 thePionPlus = G4PionPlus::PionPlus();
179 thePionMinus= G4PionMinus::PionMinus();
180
181 fTarget = G4Proton::Proton(); // later vmg
182 fMassTarg = fTarget->GetPDGMass();
183 fMassProj = fProjectile->GetPDGMass();
184 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
185 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
186 fHadrNuclXsc = new G4HadronNucleonXsc();
187
188 fEnergyBin = 200;
189 fBinT = 514; // 514; // 500; // 200;
190
191 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
192
193 fTableT = 0;
194 fOldTkin = 0.;
195
196 SetParameters();
197}
198
199
200
201//////////////////////////////////////////////////////////////////////////////
202//
203// Destructor
204
205G4hhElastic::~G4hhElastic()
206{
207 if ( fEnergyVector ) {
208 delete fEnergyVector;
209 fEnergyVector = 0;
210 }
211
212 for ( std::vector<G4PhysicsTable*>::iterator it = fBankT.begin();
213 it != fBankT.end(); ++it ) {
214 if ( (*it) ) (*it)->clearAndDestroy();
215 delete *it;
216 *it = 0;
217 }
218 fTableT = 0;
219 if(fHadrNuclXsc) delete fHadrNuclXsc;
220}
221
222/////////////////////////////////////////////////////////////////////////////
223///////////////////// Table preparation and reading ////////////////////////
224
225
226//////////////////////////////////////////////////////////////////////////////
227//
228// Initialisation for given particle on the proton target
229
230void G4hhElastic::Initialise()
231{
232 // pp,pn
233
234 fProjectile = G4Proton::Proton();
235 BuildTableT(fTarget, fProjectile);
236 fBankT.push_back(fTableT); // 0
237
238 // pi+-p
239
240 fProjectile = G4PionPlus::PionPlus();
241 BuildTableT(fTarget, fProjectile);
242 fBankT.push_back(fTableT); // 1
243 //K+-p
244 fProjectile = G4KaonPlus::KaonPlus();
245 BuildTableT(fTarget, fProjectile);
246 fBankT.push_back(fTableT); // 2
247
248}
249
250///////////////////////////////////////////////////////////////////////////////
251//
252// Build for given particle and proton table of momentum transfers.
253
254void G4hhElastic::BuildTableT( G4ParticleDefinition* target, G4ParticleDefinition* projectile) // , G4double plab)
255{
256 G4int iTkin, jTransfer;
257 G4double plab, Tkin, tMax;
258 G4double t1, t2, dt, delta = 0., sum = 0.;
259
260 fTarget = target;
261 fProjectile = projectile;
262 fMassTarg = fTarget->GetPDGMass();
263 fMassProj = fProjectile->GetPDGMass();
264 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
265 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
266
267 G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral;
268 // G4HadronNucleonXsc* hnXsc = new G4HadronNucleonXsc();
269 fTableT = new G4PhysicsTable(fEnergyBin);
270
271 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
272 {
273 Tkin = fEnergyVector->GetLowEdgeEnergy(iTkin);
274 plab = std::sqrt( Tkin*( Tkin + 2*fMassProj ) );
275 // G4DynamicParticle* theDynamicParticle = new G4DynamicParticle(projectile,
276 // G4ParticleMomentum(0.,0.,1.),
277 // Tkin);
278 // fSigmaTot = fHadrNuclXsc->GetHadronNucleonXscNS( theDynamicParticle, target );
279
280 SetParametersCMS( plab );
281
282 tMax = 4.*fPcms*fPcms;
283 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???
284
285 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);
286 sum = 0.;
287 dt = tMax/fBinT;
288
289 // for(j = 1; j < fBinT; j++)
290
291 for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer--)
292 {
293 t1 = dt*(jTransfer-1);
294 t2 = t1 + dt;
295
296 if( fMassProj > 900.*MeV ) // pp, pn
297 {
298 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
299 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
300 }
301 else // pi+-p, K+-p
302 {
303 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
304 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
305 }
306 sum += delta;
307 vectorT->PutValue( jTransfer-1, t1, sum ); // t2
308 }
309 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
310 fTableT->insertAt( iTkin, vectorT );
311 // delete theDynamicParticle;
312 }
313 // delete hnXsc;
314
315 return;
316}
317
318////////////////////////////////////////////////////////////////////////////
319//
320// Return inv momentum transfer -t > 0 from initialisation table
321
322G4double G4hhElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p,
323 G4int, G4int )
324{
325 G4int iTkin, iTransfer;
326 G4double t, t2, position, m1 = aParticle->GetPDGMass();
327 G4double Tkin = std::sqrt(m1*m1+p*p) - m1;
328
329 if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
330 {
331 fTableT = fBankT[0];
332 }
333 if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
334 {
335 fTableT = fBankT[1];
336 }
337 if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
338 {
339 fTableT = fBankT[2];
340 }
341
342 G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
343 G4double deltaMax = 1.e-2;
344
345 if ( delta < deltaMax ) iTkin = fInTkin;
346 else
347 {
348 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
349 {
350 if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
351 }
352 }
353 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
354 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
355
356 fOldTkin = Tkin;
357 fInTkin = iTkin;
358
359 if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges
360 {
361 position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
362
363 // G4cout<<"position = "<<position<<G4endl;
364
365 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
366 {
367 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
368 }
369 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
370
371 // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
372
373 t = GetTransfer(iTkin, iTransfer, position);
374
375 // G4cout<<"t = "<<t<<G4endl;
376 }
377 else // Tkin inside between energy table edges
378 {
379 // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();
380 position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
381
382 // G4cout<<"position = "<<position<<G4endl;
383
384 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
385 {
386 // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
387 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
388 }
389 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
390
391 // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
392
393 t2 = GetTransfer(iTkin, iTransfer, position);
394 return t2;
395 /*
396 G4double t1, E1, E2, W, W1, W2;
397 // G4cout<<"t2 = "<<t2<<G4endl;
398
399 E2 = fEnergyVector->GetLowEdgeEnergy(iTkin);
400
401 // G4cout<<"E2 = "<<E2<<G4endl;
402
403 iTkin--;
404
405 // position = (*(*fTableT)(iTkin))(fBinT-2)*G4UniformRand();
406
407 // G4cout<<"position = "<<position<<G4endl;
408
409 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
410 {
411 // if( position < (*(*fTableT)(iTkin))(iTransfer) ) break;
412 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
413 }
414 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
415
416 t1 = GetTransfer(iTkin, iTransfer, position);
417
418 // G4cout<<"t1 = "<<t1<<G4endl;
419
420 E1 = fEnergyVector->GetLowEdgeEnergy(iTkin);
421
422 // G4cout<<"E1 = "<<E1<<G4endl;
423
424 W = 1.0/(E2 - E1);
425 W1 = (E2 - Tkin)*W;
426 W2 = (Tkin - E1)*W;
427
428 t = W1*t1 + W2*t2;
429 */
430 }
431 return t;
432}
433
434
435////////////////////////////////////////////////////////////////////////////
436//
437// Return inv momentum transfer -t > 0 from initialisation table
438
439G4double G4hhElastic::SampleBisectionalT( const G4ParticleDefinition* aParticle, G4double p)
440{
441 G4int iTkin, iTransfer;
442 G4double t, position, m1 = aParticle->GetPDGMass();
443 G4double Tkin = std::sqrt(m1*m1+p*p) - m1;
444
445 if( aParticle == G4Proton::Proton() || aParticle == G4Neutron::Neutron() )
446 {
447 fTableT = fBankT[0];
448 }
449 if( aParticle == G4PionPlus::PionPlus() || aParticle == G4PionMinus::PionMinus() )
450 {
451 fTableT = fBankT[1];
452 }
453 if( aParticle == G4KaonPlus::KaonPlus() || aParticle == G4KaonMinus::KaonMinus() )
454 {
455 fTableT = fBankT[2];
456 }
457 G4double delta = std::abs(Tkin - fOldTkin)/(Tkin + fOldTkin);
458 G4double deltaMax = 1.e-2;
459
460 if ( delta < deltaMax ) iTkin = fInTkin;
461 else
462 {
463 for( iTkin = 0; iTkin < fEnergyBin; iTkin++ )
464 {
465 if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
466 }
467 }
468 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
469 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
470
471 fOldTkin = Tkin;
472 fInTkin = iTkin;
473
474 if (iTkin == fEnergyBin -1 || iTkin == 0 ) // the table edges
475 {
476 position = (*(*fTableT)(iTkin))(0)*G4UniformRand();
477
478 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
479 {
480 if( position >= (*(*fTableT)(iTkin))(iTransfer) ) break;
481 }
482 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
483
484 t = GetTransfer(iTkin, iTransfer, position);
485
486
487 }
488 else // Tkin inside between energy table edges
489 {
490 G4double rand = G4UniformRand();
491 position = (*(*fTableT)(iTkin))(0)*rand;
492
493 //
494 // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBinT-2);
495 G4int sTransfer = 0, fTransfer = fBinT - 2, dTransfer = fTransfer - sTransfer;
496 G4double y2;
497
498 for( iTransfer = 0; iTransfer < fBinT - 1; iTransfer++ )
499 {
500 // dTransfer %= 2;
501 dTransfer /= 2;
502 // dTransfer *= 0.5;
503 y2 = (*(*fTableT)(iTkin))( sTransfer + dTransfer );
504
505 if( y2 > position ) sTransfer += dTransfer;
506
507 // if( dTransfer <= 1 ) break;
508 if( dTransfer < 1 ) break;
509 }
510 t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sTransfer); // +(-0.5+rand)*(*fTableT)(iTkin)->GetLowEdgeEnergy(3);
511 }
512 return t;
513}
514
515
516///////////////////////////////////////////////////////////////////////////////
517//
518// Build for given particle and proton table of momentum transfers.
519
520void G4hhElastic::BuildTableTest( G4ParticleDefinition* target, G4ParticleDefinition* projectile, G4double plab)
521{
522 G4int jTransfer;
523 G4double tMax; // , sQq, sQG;
524 G4double t1, t2, dt, delta = 0., sum = 0. ; // , threshold;
525
526 fTarget = target;
527 fProjectile = projectile;
528 fMassTarg = fTarget->GetPDGMass();
529 fMassProj = fProjectile->GetPDGMass();
530 fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg+fMassProj);
531 fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg-fMassProj);
532 fSpp = fMassProj*fMassProj + fMassTarg*fMassTarg + 2.*fMassTarg*std::sqrt(plab*plab + fMassProj*fMassProj);
533 fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp - fMassDif2)/4./fSpp);
534
535 G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fMassProj = "<<fMassProj<<" MeV"<<G4endl;
536 tMax = 4.*fPcms*fPcms;
537 if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV; // Check vs. energy ???
538
539
540 G4Integrator<G4hhElastic,G4double(G4hhElastic::*)(G4double)> integral;
541 fTableT = new G4PhysicsTable(1);
542 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fBinT-1);
543
544 sum = 0.;
545 dt = tMax/G4double(fBinT);
546 G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV; fPcms = "<<fPcms/GeV
547 <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt = "<<dt/GeV/GeV<<" GeV2"<<G4endl;
548
549 // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fRB*GeV<<G4endl;
550
551 // for(jTransfer = 1; jTransfer < fBinT; jTransfer++)
552 for( jTransfer = fBinT-1; jTransfer >= 1; jTransfer-- )
553 {
554 t1 = dt*(jTransfer-1);
555 t2 = t1 + dt;
556
557 if( fMassProj > 900.*MeV ) // pp, pn
558 {
559 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123, t1, t2);
560 // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, tMax);
561 }
562 else // pi+-p, K+-p
563 {
564 delta = integral.Legendre10(this, &G4hhElastic::GetdsdtF123qQgG, t1, t2);
565 // threshold = integral.Legendre96(this, &G4hhElastic::GetdsdtF123qQgG, t1, tMax);
566 // delta = integral.Legendre96(this, &G4hhElastic::GetdsdtF123, t1, t2);
567 }
568 sum += delta;
569 // G4cout<<delta<<"\t"<<sum<<"\t"<<threshold<<G4endl;
570
571 // sQq = GetdsdtF123(q1);
572 // sQG = GetdsdtF123qQgG(q1);
573 // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/millibarn<<"\t"<<sQq*GeV*GeV/millibarn<<G4endl;
574 // G4cout<<"sum = "<<sum<<", ";
575
576 vectorT->PutValue( jTransfer-1, t1, sum ); // t2
577 }
578 // vectorT->PutValue( fBinT-1, dt*(fBinT-1), 0. ); // t2
579 fTableT->insertAt( 0, vectorT );
580 fBankT.push_back( fTableT ); // 0
581
582 // for(jTransfer = 0; jTransfer < fBinT-1; jTransfer++)
583 // G4cout<<(*(*fTableT)(0))(jTransfer)/sum<<"\t\t"<<G4Pow::GetInstance()->powN(2.,-jTransfer)<<G4endl;
584
585 return;
586}
587
588
589////////////////////////////////////////////////////////////////////////////
590//
591// Return inv momentum transfer -t > 0 from initialisation table
592
593G4double G4hhElastic::SampleTest(G4double tMin ) // const G4ParticleDefinition* aParticle, )
594{
595 G4int iTkin, iTransfer, iTmin;
596 G4double t, position;
597 // G4double qMin = std::sqrt(tMin);
598
599 fTableT = fBankT[0];
600 iTkin = 0;
601
602 for(iTransfer = 0; iTransfer < fBinT-1; iTransfer++)
603 {
604 // if( qMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
605 if( tMin <= (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer) ) break;
606 }
607 iTmin = iTransfer-1;
608 if(iTmin < 0 ) iTmin = 0;
609
610 position = (*(*fTableT)(iTkin))(iTmin)*G4UniformRand();
611
612 for( iTmin = 0; iTransfer < fBinT-1; iTransfer++)
613 {
614 if( position > (*(*fTableT)(iTkin))(iTransfer) ) break;
615 }
616 if (iTransfer >= fBinT-1) iTransfer = fBinT-2;
617
618 t = GetTransfer(iTkin, iTransfer, position);
619
620 return t;
621}
622
623
624/////////////////////////////////////////////////////////////////////////////////
625//
626// Check with PAI sampling
627
628G4double
629G4hhElastic:: GetTransfer( G4int iTkin, G4int iTransfer, G4double position )
630{
631 G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;
632
633 if( iTransfer == 0 )
634 {
635 randTransfer = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);
636 // iTransfer++;
637 }
638 else
639 {
640 if ( iTransfer >= G4int((*fTableT)(iTkin)->GetVectorLength()) )
641 {
642 iTransfer = G4int((*fTableT)(iTkin)->GetVectorLength() - 1);
643 }
644 y1 = (*(*fTableT)(iTkin))(iTransfer-1);
645 y2 = (*(*fTableT)(iTkin))(iTransfer);
646
647 x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer-1);
648 x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(iTransfer);
649
650 delta = y2 - y1;
651 mean = y2 + y1;
652
653 if ( x1 == x2 ) randTransfer = x2;
654 else
655 {
656 // if ( y1 == y2 )
657 if ( delta < epsilon*mean )
658 randTransfer = x1 + ( x2 - x1 )*G4UniformRand();
659 else randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
660 }
661 }
662 return randTransfer;
663}
664
665const G4double G4hhElastic::theNuclNuclData[19][6] =
666{
667 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1/GeV, fBq, fBQ, fImCof
668
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 }, // pp 4GeV/c
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 }, // np 9 GeV/c
673 { 4.62126, 6, 6, 0.03, 0.770111, 6.5 }, // pp 10.4 GeV/c
674
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 }, // pp 19.2 GeV/c
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 }, // pp 50 GeV/c
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 }, // np 57.5 GeV/c
681
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.5 }, // pp 200 GeV/c
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, 0.576215, 7.5 }, // np 309.5 GeV/c
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 }, // pp 58.2 GeV
688
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 }, // pb-p 1960 GeV
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 } // pp TOTEM 13 TeV
693
694};
695
696//////////////////////////////////////////////////////////////////////////////////
697
698const G4double G4hhElastic::thePiKaNuclData[8][6] =
699{
700 // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1/GeV, fBq, fBQ, fImCof
701
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.250601, 1.3 }, // pipp 4.122 GeV/c
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 }, // pipp 32 GeV/c
706 { 9.7328, 5, 4.5, 0.05, 0.416319, 5.5 }, // pipp 50 GeV/c
707
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 }, // pipp 147 GeV/c
710 { 19.3961, 5, 4.5, 0.05, 0.413477, 3.5 } // pimp 200 GeV/c
711
712};
713
714//
715//
716/////////////////////////////////////////////////////////////////////////////////
epsilon
G4double epsilon(G4double density, G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:95
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4HadronicInteraction::SetMinEnergy
void SetMinEnergy(G4double anEnergy)
Definition: G4HadronicInteraction.hh:89
G4HadronicInteraction::SetMaxEnergy
void SetMaxEnergy(const G4double anEnergy)
Definition: G4HadronicInteraction.hh:102
G4HadronicParameters::Instance
static G4HadronicParameters * Instance()
Definition: G4HadronicParameters.cc:54
G4KaonMinus::KaonMinus
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:112
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:112
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
G4PhysicsFreeVector::PutValue
void PutValue(const std::size_t index, const G4double e, const G4double value)
Definition: G4PhysicsFreeVector.hh:95
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(const std::size_t index) const
G4PionMinus::PionMinus
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:97
G4PionPlus::PionPlus
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:97
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:92
G4hhElastic::GetdsdtF123
G4double GetdsdtF123(G4double q)
Definition: G4hhElastic.hh:1128
G4hhElastic::SampleBisectionalT
G4double SampleBisectionalT(const G4ParticleDefinition *p, G4double plab)
Definition: G4hhElastic.cc:439
G4hhElastic::BuildTableTest
void BuildTableTest(G4ParticleDefinition *target, G4ParticleDefinition *projectile, G4double plab)
Definition: G4hhElastic.cc:520
G4hhElastic::SampleInvariantT
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int, G4int)
Definition: G4hhElastic.cc:322
G4hhElastic::SetParametersCMS
void SetParametersCMS(G4double plab)
Definition: G4hhElastic.hh:316
G4hhElastic::Initialise
void Initialise()
Definition: G4hhElastic.cc:230
G4hhElastic::GetdsdtF123qQgG
G4double GetdsdtF123qQgG(G4double q)
Definition: G4hhElastic.hh:746
G4hhElastic::SampleTest
G4double SampleTest(G4double tMin)
Definition: G4hhElastic.cc:593
G4hhElastic::~G4hhElastic
virtual ~G4hhElastic()
Definition: G4hhElastic.cc:205
G4hhElastic::BuildTableT
void BuildTableT(G4ParticleDefinition *target, G4ParticleDefinition *projectile)
Definition: G4hhElastic.cc:254
G4hhElastic::GetTransfer
G4double GetTransfer(G4int iMomentum, G4int iTransfer, G4double position)
Definition: G4hhElastic.cc:629
G4hhElastic::SetParameters
void SetParameters()
Definition: G4hhElastic.hh:271