Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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 // $Id: G4hhElastic.cc,v 1.5 2010-11-09 09:04:29 grichine Exp $
27 // GEANT4 tag $Name: not supported by cvs2svn $
28 //
29 //
30 // Physics model class G4hhElastic
31 //
32 //
33 // G4 Model: qQ hadron hadron elastic scattering with 4-momentum balance
34 //
35 // 02.05.2014 V. Grichine 1-st version
36 //
37 
38 #include "G4hhElastic.hh"
39 #include "G4ParticleTable.hh"
40 #include "G4ParticleDefinition.hh"
41 #include "G4IonTable.hh"
42 #include "G4NucleiProperties.hh"
43 
44 #include "Randomize.hh"
45 #include "G4Integrator.hh"
46 #include "globals.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 
50 #include "G4Proton.hh"
51 #include "G4Neutron.hh"
52 #include "G4PionPlus.hh"
53 #include "G4PionMinus.hh"
54 
55 #include "G4Element.hh"
56 #include "G4ElementTable.hh"
57 #include "G4PhysicsTable.hh"
58 #include "G4PhysicsLogVector.hh"
59 #include "G4PhysicsFreeVector.hh"
60 
61 #include "G4HadronNucleonXsc.hh"
62 
63 #include "G4Pow.hh"
64 
65 using namespace std;
66 
67 
69 //
70 // Tracking constructor. Target is proton
71 
72 
74  : G4HadronElastic("HadrHadrElastic")
75 {
76  SetMinEnergy( 1.*GeV );
77  SetMaxEnergy( 10000.*TeV );
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 
110 //
111 // test constructor
112 
113 
115  : G4HadronElastic("HadrHadrElastic")
116 {
117  SetMinEnergy( 1.*GeV );
118  SetMaxEnergy( 10000.*TeV );
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 
155 //
156 // constructor used for low mass diffraction
157 
158 
160  : G4HadronElastic("HadrHadrElastic")
161 {
162  SetMinEnergy( 1.*GeV );
163  SetMaxEnergy( 10000.*TeV );
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 
202 //
203 // Destructor
204 
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 
224 
225 
227 //
228 // Initialisation for given particle on the proton target
229 
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 
251 //
252 // Build for given particle and proton table of momentum transfers.
253 
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 
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 
319 //
320 // Return inv momentum transfer -t > 0 from initialisation table
321 
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 
436 //
437 // Return inv momentum transfer -t > 0 from initialisation table
438 
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 
517 //
518 // Build for given particle and proton table of momentum transfers.
519 
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 
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 
590 //
591 // Return inv momentum transfer -t > 0 from initialisation table
592 
593 G4double 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 
625 //
626 // Check with PAI sampling
627 
628 G4double
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 = (*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 
665 const G4double G4hhElastic::theNuclNuclData[18][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
692 
693 };
694 
696 
697 const G4double G4hhElastic::thePiKaNuclData[8][6] =
698 {
699  // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1/GeV, fBq, fBQ, fImCof
700 
701  { 2.5627, 3.8, 3.3, 0.22, 0.222, 1.5 }, // pipp 3.017 GeV/c
702  { 2.93928, 4.3, 3.8, 0.2, 0.250601, 1.3 }, // pipp 4.122 GeV/c
703  { 3.22326, 4.8, 4.3, 0.13, 0.32751, 2.5 }, // pipp 5.055 GeV/c
704  { 7.80704, 5.5, 5, 0.13, 0.340631, 2.5 }, // pipp 32 GeV/c
705  { 9.7328, 5, 4.5, 0.05, 0.416319, 5.5 }, // pipp 50 GeV/c
706 
707  { 13.7315, 5.3, 4.8, 0.05, 0.418426, 5.5 }, // pipp 100 GeV/c
708  { 16.6359, 6.3, 5.8, 0.05, 0.423817, 5.5 }, // pipp 147 GeV/c
709  { 19.3961, 5, 4.5, 0.05, 0.413477, 3.5 } // pimp 200 GeV/c
710 
711 };
712 
713 //
714 //
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
void SetParameters()
Definition: G4hhElastic.hh:273
void SetParametersCMS(G4double plab)
Definition: G4hhElastic.hh:318
virtual ~G4hhElastic()
Definition: G4hhElastic.cc:205
const XML_Char * target
Definition: expat.h:268
const char * p
Definition: xmltok.h:285
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
static constexpr double TeV
Definition: G4SIunits.hh:218
void BuildTableT(G4ParticleDefinition *target, G4ParticleDefinition *projectile)
Definition: G4hhElastic.cc:254
void SetMinEnergy(G4double anEnergy)
void BuildTableTest(G4ParticleDefinition *target, G4ParticleDefinition *projectile, G4double plab)
Definition: G4hhElastic.cc:520
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double GetTransfer(G4int iMomentum, G4int iTransfer, G4double position)
Definition: G4hhElastic.cc:629
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4double GetdsdtF123qQgG(G4double q)
Definition: G4hhElastic.hh:748
G4double GetdsdtF123(G4double q)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int, G4int)
Definition: G4hhElastic.cc:322
G4double GetPDGMass() const
G4double SampleBisectionalT(const G4ParticleDefinition *p, G4double plab)
Definition: G4hhElastic.cc:439
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
tuple t1
Definition: plottest35.py:33
static constexpr double GeV
Definition: G4SIunits.hh:217
void insertAt(size_t, G4PhysicsVector *)
void SetMaxEnergy(const G4double anEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
void Initialise()
Definition: G4hhElastic.cc:230
double G4double
Definition: G4Types.hh:76
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
static constexpr double keV
Definition: G4SIunits.hh:216
double epsilon(double density, double temperature)
G4double SampleTest(G4double tMin)
Definition: G4hhElastic.cc:593