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 ]

Diff markup

Differences between /processes/hadronic/models/coherent_elastic/src/G4hhElastic.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4hhElastic.cc (Version 9.6.p2)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 //                                                
 28 // Physics model class G4hhElastic                
 29 //                                                
 30 //                                                
 31 // G4 Model: qQ hadron hadron elastic scatteri    
 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    
 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( the    
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    
115   : G4HadronElastic("HadrHadrElastic")            
116 {                                                 
117   SetMinEnergy( 1.*GeV );                         
118   SetMaxEnergy( G4HadronicParameters::Instance    
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)*(fMassTa    
138   fMassDif2   = (fMassTarg-fMassProj)*(fMassTa    
139   fHadrNuclXsc = new G4HadronNucleonXsc();        
140                                                   
141   fEnergyBin = 200;                               
142   fBinT      = 514; // 200;                       
143                                                   
144   fEnergyVector =  new G4PhysicsLogVector( the    
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    
160   : G4HadronElastic("HadrHadrElastic")            
161 {                                                 
162   SetMinEnergy( 1.*GeV );                         
163   SetMaxEnergy( G4HadronicParameters::Instance    
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)*(fMassTa    
185   fMassDif2   = (fMassTarg-fMassProj)*(fMassTa    
186   fHadrNuclXsc = new G4HadronNucleonXsc();        
187                                                   
188   fEnergyBin = 200;                               
189   fBinT  = 514; // 514; // 500; // 200;           
190                                                   
191   fEnergyVector =  new G4PhysicsLogVector( the    
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    
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 r    
224                                                   
225                                                   
226 //////////////////////////////////////////////    
227 //                                                
228 // Initialisation for given particle on the pr    
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 o    
253                                                   
254 void G4hhElastic::BuildTableT( G4ParticleDefin    
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)*(fMassTa    
265   fMassDif2   = (fMassTarg-fMassProj)*(fMassTa    
266                                                   
267   G4Integrator<G4hhElastic,G4double(G4hhElasti    
268   // G4HadronNucleonXsc*          hnXsc = new     
269   fTableT = new G4PhysicsTable(fEnergyBin);       
270                                                   
271   for( iTkin = 0; iTkin < fEnergyBin; iTkin++)    
272   {                                               
273     Tkin  = fEnergyVector->GetLowEdgeEnergy(iT    
274     plab  = std::sqrt( Tkin*( Tkin + 2*fMassPr    
275     // G4DynamicParticle*  theDynamicParticle     
276     //                                            
277     //                                            
278     // fSigmaTot = fHadrNuclXsc->GetHadronNucl    
279                                                   
280     SetParametersCMS( plab );                     
281                                                   
282     tMax  = 4.*fPcms*fPcms;                       
283     if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*Ge    
284                                                   
285     G4PhysicsFreeVector* vectorT = new G4Physi    
286     sum = 0.;                                     
287     dt  = tMax/fBinT;                             
288                                                   
289     // for(j = 1; j < fBinT; j++)                 
290                                                   
291     for( jTransfer = fBinT-1; jTransfer >= 1;     
292     {                                             
293       t1 = dt*(jTransfer-1);                      
294       t2 = t1 + dt;                               
295                                                   
296       if( fMassProj > 900.*MeV ) // pp, pn        
297       {                                           
298         delta = integral.Legendre10(this, &G4h    
299         // delta = integral.Legendre96(this, &    
300       }                                           
301       else // pi+-p, K+-p                         
302       {                                           
303         delta = integral.Legendre10(this, &G4h    
304         // delta = integral.Legendre96(this, &    
305       }                                           
306       sum += delta;                               
307       vectorT->PutValue( jTransfer-1, t1, sum     
308     }                                             
309     // vectorT->PutValue( fBinT-1, dt*(fBinT-1    
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 in    
321                                                   
322 G4double G4hhElastic::SampleInvariantT( const     
323                                                   
324 {                                                 
325   G4int iTkin, iTransfer;                         
326   G4double t, t2, position, m1 = aParticle->Ge    
327   G4double Tkin = std::sqrt(m1*m1+p*p) - m1;      
328                                                   
329   if( aParticle == G4Proton::Proton() || aPart    
330   {                                               
331     fTableT = fBankT[0];                          
332   }                                               
333   if( aParticle == G4PionPlus::PionPlus() || a    
334   {                                               
335     fTableT = fBankT[1];                          
336   }                                               
337   if( aParticle == G4KaonPlus::KaonPlus() || a    
338   {                                               
339     fTableT = fBankT[2];                          
340   }                                               
341                                                   
342   G4double delta    = std::abs(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->GetLowEdgeEner    
351     }                                             
352   }                                               
353   if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi    
354   if ( iTkin < 0 )           iTkin = 0; // aga    
355                                                   
356   fOldTkin = Tkin;                                
357   fInTkin = iTkin;                                
358                                                   
359   if (iTkin == fEnergyBin -1 || iTkin == 0 )      
360   {                                               
361     position = (*(*fTableT)(iTkin))(0)*G4Unifo    
362                                                   
363     // G4cout<<"position = "<<position<<G4endl    
364                                                   
365     for(iTransfer = 0; iTransfer < fBinT-1; iT    
366     {                                             
367       if( position >= (*(*fTableT)(iTkin))(iTr    
368     }                                             
369     if (iTransfer >= fBinT-1) iTransfer = fBin    
370                                                   
371     // G4cout<<"iTransfer = "<<iTransfer<<G4en    
372                                                   
373     t = GetTransfer(iTkin, iTransfer, position    
374                                                   
375     // G4cout<<"t = "<<t<<G4endl;                 
376   }                                               
377   else  // Tkin inside between energy table ed    
378   {                                               
379     // position = (*(*fTableT)(iTkin))(fBinT-2    
380     position = (*(*fTableT)(iTkin))(0)*G4Unifo    
381                                                   
382     // G4cout<<"position = "<<position<<G4endl    
383                                                   
384     for(iTransfer = 0; iTransfer < fBinT-1; iT    
385     {                                             
386       // if( position < (*(*fTableT)(iTkin))(i    
387       if( position >= (*(*fTableT)(iTkin))(iTr    
388     }                                             
389     if (iTransfer >= fBinT-1) iTransfer = fBin    
390                                                   
391     // G4cout<<"iTransfer = "<<iTransfer<<G4en    
392                                                   
393     t2  = GetTransfer(iTkin, iTransfer, positi    
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    
406                                                   
407     // G4cout<<"position = "<<position<<G4endl    
408                                                   
409     for(iTransfer = 0; iTransfer < fBinT-1; iT    
410     {                                             
411       // if( position < (*(*fTableT)(iTkin))(i    
412       if( position >= (*(*fTableT)(iTkin))(iTr    
413     }                                             
414     if (iTransfer >= fBinT-1) iTransfer = fBin    
415                                                   
416     t1  = GetTransfer(iTkin, iTransfer, positi    
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 in    
438                                                   
439 G4double G4hhElastic::SampleBisectionalT( cons    
440 {                                                 
441   G4int iTkin, iTransfer;                         
442   G4double t,  position, m1 = aParticle->GetPD    
443   G4double Tkin = std::sqrt(m1*m1+p*p) - m1;      
444                                                   
445   if( aParticle == G4Proton::Proton() || aPart    
446   {                                               
447     fTableT = fBankT[0];                          
448   }                                               
449   if( aParticle == G4PionPlus::PionPlus() || a    
450   {                                               
451     fTableT = fBankT[1];                          
452   }                                               
453   if( aParticle == G4KaonPlus::KaonPlus() || a    
454   {                                               
455     fTableT = fBankT[2];                          
456   }                                               
457   G4double delta    = std::abs(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->GetLowEdgeEner    
466     }                                             
467   }                                               
468   if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi    
469   if ( iTkin < 0 )           iTkin = 0; // aga    
470                                                   
471   fOldTkin = Tkin;                                
472   fInTkin = iTkin;                                
473                                                   
474   if (iTkin == fEnergyBin -1 || iTkin == 0 )      
475   {                                               
476     position = (*(*fTableT)(iTkin))(0)*G4Unifo    
477                                                   
478     for(iTransfer = 0; iTransfer < fBinT-1; iT    
479     {                                             
480       if( position >= (*(*fTableT)(iTkin))(iTr    
481     }                                             
482     if (iTransfer >= fBinT-1) iTransfer = fBin    
483                                                   
484     t = GetTransfer(iTkin, iTransfer, position    
485                                                   
486                                                   
487   }                                               
488   else  // Tkin inside between energy table ed    
489   {                                               
490     G4double rand = G4UniformRand();              
491     position = (*(*fTableT)(iTkin))(0)*rand;      
492                                                   
493     //                                            
494     // (*fTableT)(iTkin)->GetLowEdgeEnergy(fBi    
495     G4int sTransfer = 0, fTransfer =  fBinT -     
496     G4double y2;                                  
497                                                   
498     for( iTransfer = 0; iTransfer < fBinT - 1;    
499     {                                             
500       // dTransfer %= 2;                          
501       dTransfer /= 2;                             
502       // dTransfer *= 0.5;                        
503       y2 = (*(*fTableT)(iTkin))( sTransfer + d    
504                                                   
505       if( y2 > position ) sTransfer += dTransf    
506                                                   
507       // if( dTransfer <= 1 ) break;              
508       if( dTransfer < 1 ) break;                  
509     }                                             
510     t = (*fTableT)(iTkin)->GetLowEdgeEnergy(sT    
511   }                                               
512   return t;                                       
513 }                                                 
514                                                   
515                                                   
516 //////////////////////////////////////////////    
517 //                                                
518 // Build for given particle and proton table o    
519                                                   
520 void G4hhElastic::BuildTableTest( G4ParticleDe    
521 {                                                 
522   G4int jTransfer;                                
523   G4double tMax; // , sQq, sQG;                   
524   G4double t1, t2, dt, delta = 0., sum = 0. ;     
525                                                   
526   fTarget     = target;                           
527   fProjectile = projectile;                       
528   fMassTarg =  fTarget->GetPDGMass();             
529   fMassProj =  fProjectile->GetPDGMass();         
530   fMassSum2 = (fMassTarg+fMassProj)*(fMassTarg    
531   fMassDif2 = (fMassTarg-fMassProj)*(fMassTarg    
532   fSpp  = fMassProj*fMassProj + fMassTarg*fMas    
533   fPcms = std::sqrt( (fSpp - fMassSum2)*(fSpp     
534                                                   
535   G4cout<<"fMassTarg = "<<fMassTarg<<" MeV; fM    
536   tMax = 4.*fPcms*fPcms;                          
537   if( tMax > 15.*GeV*GeV ) tMax = 15.*GeV*GeV;    
538                                                   
539                                                   
540   G4Integrator<G4hhElastic,G4double(G4hhElasti    
541   fTableT = new G4PhysicsTable(1);                
542   G4PhysicsFreeVector* vectorT = new G4Physics    
543                                                   
544   sum = 0.;                                       
545   dt = tMax/G4double(fBinT);                      
546   G4cout<<"s = "<<std::sqrt(fSpp)/GeV<<" GeV;     
547   <<" GeV; qMax = "<<tMax/GeV/GeV<<" GeV2; dt     
548                                                   
549   // G4cout<<"fRA = "<<fRA*GeV<<"; fRB = "<<fR    
550                                                   
551     // for(jTransfer = 1; jTransfer < fBinT; j    
552   for( jTransfer = fBinT-1; jTransfer >= 1; jT    
553   {                                               
554     t1 = dt*(jTransfer-1);                        
555     t2 = t1 + dt;                                 
556                                                   
557     if( fMassProj > 900.*MeV ) // pp, pn          
558     {                                             
559       delta = integral.Legendre10(this, &G4hhE    
560       // threshold = integral.Legendre96(this,    
561     }                                             
562     else // pi+-p, K+-p                           
563     {                                             
564       delta = integral.Legendre10(this, &G4hhE    
565       // threshold = integral.Legendre96(this,    
566       // delta = integral.Legendre96(this, &G4    
567     }                                             
568     sum += delta;                                 
569     // G4cout<<delta<<"\t"<<sum<<"\t"<<thresho    
570                                                   
571     // sQq = GetdsdtF123(q1);                     
572     // sQG = GetdsdtF123qQgG(q1);                 
573     // G4cout<<q1/GeV<<"\t"<<sQG*GeV*GeV/milli    
574     // G4cout<<"sum = "<<sum<<", ";               
575                                                   
576     vectorT->PutValue( jTransfer-1, t1, sum );    
577   }                                               
578   // vectorT->PutValue( fBinT-1, dt*(fBinT-1),    
579   fTableT->insertAt( 0, vectorT );                
580   fBankT.push_back( fTableT );  // 0              
581                                                   
582   // for(jTransfer = 0; jTransfer < fBinT-1; j    
583   //   G4cout<<(*(*fTableT)(0))(jTransfer)/sum    
584                                                   
585   return;                                         
586 }                                                 
587                                                   
588                                                   
589 //////////////////////////////////////////////    
590 //                                                
591 // Return inv momentum transfer -t > 0 from in    
592                                                   
593 G4double G4hhElastic::SampleTest(G4double tMin    
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; iTra    
603   {                                               
604     // if( qMin <= (*fTableT)(iTkin)->GetLowEd    
605     if( tMin <= (*fTableT)(iTkin)->GetLowEdgeE    
606   }                                               
607   iTmin = iTransfer-1;                            
608   if(iTmin < 0 ) iTmin = 0;                       
609                                                   
610   position = (*(*fTableT)(iTkin))(iTmin)*G4Uni    
611                                                   
612   for( iTmin = 0; iTransfer < fBinT-1; iTransf    
613   {                                               
614     if( position > (*(*fTableT)(iTkin))(iTrans    
615   }                                               
616   if (iTransfer >= fBinT-1) iTransfer = fBinT-    
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     
630 {                                                 
631   G4double x1, x2, y1, y2, randTransfer, delta    
632                                                   
633   if( iTransfer == 0 )                            
634   {                                               
635     randTransfer = (*fTableT)(iTkin)->GetLowEd    
636     // iTransfer++;                               
637   }                                               
638   else                                            
639   {                                               
640     if ( iTransfer >= G4int((*fTableT)(iTkin)-    
641     {                                             
642       iTransfer = G4int((*fTableT)(iTkin)->Get    
643     }                                             
644     y1 = (*(*fTableT)(iTkin))(iTransfer-1);       
645     y2 = (*(*fTableT)(iTkin))(iTransfer);         
646                                                   
647     x1 = (*fTableT)(iTkin)->GetLowEdgeEnergy(i    
648     x2 = (*fTableT)(iTkin)->GetLowEdgeEnergy(i    
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 )*G4U    
659       else randTransfer = x1 + ( position - y1    
660     }                                             
661   }                                               
662   return randTransfer;                            
663 }                                                 
664                                                   
665 const G4double G4hhElastic::theNuclNuclData[19    
666 {                                                 
667   // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1    
668                                                   
669   { 2.76754,  4.8,  4.8,  0.05, 0.742441, 10.5    
670   { 3.07744,  5.4,  5.4,  0.02, 0.83818,  6.5     
671   { 3.36305,  5.2,  5.2,  0.02, 0.838893, 7.5     
672   { 4.32941,  6,  6,  0.03, 0.769389, 7.5  },     
673   { 4.62126,  6,  6,  0.03, 0.770111, 6.5  },     
674                                                   
675   { 5.47416,  4.5,  4.5,  0.03, 0.813185, 7.5     
676   { 6.15088,  6.5,  6.5,  0.02, 0.799539, 6.5     
677   { 6.77474,  5.2,  5.2,  0.03, 0.784901, 7.5     
678   { 9.77775,  7,  7,  0.03, 0.742531, 6.5  },     
679                 // {9.77775,  7,  7,  0.011,      
680   { 10.4728,  5.2,  5.2,  0.03, 0.780439, 7.5     
681                                                   
682   { 13.7631,  7,  7,  0.008,  0.8664,             
683   { 19.4184,  6.8,  6.8,  0.009,  0.861337, 2.    
684   { 23.5, 6.8,  6.8,  0.007,  0.878112, 1.5  }    
685                 // {24.1362,  6.4,  6.4,  0.09    
686   { 24.1362,  7.2,  7.2,  0.008,  0.864745, 5.    
687   { 52.8, 6.8,  6.8,  0.008,  0.871929, 1.5  }    
688                                                   
689   { 546,        7.4,  7.4,  0.013,  0.845877,     
690   { 1960, 7.8,  7.8,  0.022,  0.809062, 7.5  }    
691   { 7000, 8,  8,  0.024,  0.820441, 5.5  },  /    
692   { 13000,  8.5,  8.5,  0.03, 0.796721, 10.5      
693                                                   
694 };                                                
695                                                   
696 //////////////////////////////////////////////    
697                                                   
698 const G4double G4hhElastic::thePiKaNuclData[8]    
699 {                                                 
700   // sqrt(fSpp) in GeV, fRA in 1/GeV, fRB in 1    
701                                                   
702   { 2.5627,     3.8,    3.3,    0.22,   0.222,    
703   { 2.93928,    4.3,    3.8,    0.2,    0.2506    
704   { 3.22326,  4.8,  4.3,  0.13, 0.32751,  2.5     
705   { 7.80704,  5.5,  5,  0.13, 0.340631, 2.5 },    
706   { 9.7328, 5,  4.5,  0.05, 0.416319, 5.5  },     
707                                                   
708   { 13.7315,  5.3,  4.8,  0.05, 0.418426, 5.5     
709   { 16.6359,  6.3,  5.8,  0.05, 0.423817, 5.5     
710   { 19.3961,  5,  4.5,  0.05, 0.413477, 3.5  }    
711                                                   
712 };                                                
713                                                   
714 //                                                
715 //                                                
716 //////////////////////////////////////////////    
717