Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4hhElastic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 65 using namespace std;
 66 
 67 
 68 /////////////////////////////////////////////////////////////////////////
 69 //
 70 // Tracking constructor. Target is proton
 71 
 72 
 73 G4hhElastic::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 
114 G4hhElastic::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 
159 G4hhElastic::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 
205 G4hhElastic::~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 
230 void 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 
254 void 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 
322 G4double 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 
439 G4double 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 
520 void 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 
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 
624 /////////////////////////////////////////////////////////////////////////////////
625 //
626 // Check with PAI sampling
627 
628 G4double 
629 G4hhElastic:: 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 
665 const 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 
698 const 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 /////////////////////////////////////////////////////////////////////////////////
717