source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4ElasticHadrNucleusHE.cc @ 962

Last change on this file since 962 was 962, checked in by garnier, 15 years ago

update processes
File size: 42.6 KB
Line 
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// $Id: G4ElasticHadrNucleusHE.cc,v 1.79 2008/01/14 10:39:13 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31//  The generator of high energy hadron-nucleus elastic scattering
32//  The hadron kinetic energy T > 1 GeV
33//  N.  Starkov 2003.
34//
35//  19.05.04 Variant for G4 6.1: The 'ApplyYourself' was changed
36//  19.11.05 The HE elastic scattering on proton is added (N.Starkov)
37//  16.11.06 The low energy boundary is shifted to T = 400 MeV (N.Starkov)
38//  23.11.06 General cleanup, ONQ0=3, use pointer instead of particle name (VI)
39//  02.05.07 Scale sampled t as p^2 (VI)
40//  15.05.07 Redesign and cleanup (V.Ivanchenko)
41//  17.05.07 cleanup (V.Grichine)
42//
43
44#include  "G4ElasticHadrNucleusHE.hh"
45#include  "Randomize.hh"
46#include  "G4ios.hh"
47#include  "G4ParticleTable.hh"
48#include  "G4IonTable.hh"
49#include  "G4Proton.hh"
50#include  "G4NistManager.hh"
51
52using namespace std;
53
54///////////////////////////////////////////////////////////////
55//
56//
57
58G4ElasticData::G4ElasticData(const G4ParticleDefinition* p, 
59                             G4int Z, G4double AWeight, G4double* eGeV)
60{ 
61  hadr     = p;
62  massGeV  = p->GetPDGMass()/GeV;
63  mass2GeV2= massGeV*massGeV;
64  AtomicWeight = G4int(AWeight + 0.5);
65
66  DefineNucleusParameters(AWeight);
67
68  limitQ2 = 35./(R1*R1);     //  (GeV/c)^2
69
70  G4double dQ2 = limitQ2/(ONQ2 - 1.);
71
72  TableQ2[0] = 0.0;
73
74  for(G4int ii = 1; ii < ONQ2; ii++) 
75  {
76    TableQ2[ii] = TableQ2[ii-1]+dQ2;
77  }
78
79  massA  = AWeight*amu_c2/GeV;
80  massA2 = massA*massA; 
81
82  for(G4int kk = 0; kk < NENERGY; kk++) 
83  {
84    dnkE[kk] = 0;
85    G4double elab = eGeV[kk] + massGeV;
86    G4double plab2= eGeV[kk]*(eGeV[kk] + 2.0*massGeV);
87    G4double Q2m  = 4.0*plab2*massA2/(mass2GeV2 + massA2 + 2.*massA2*elab);
88
89    if(Z == 1 && p == G4Proton::Proton()) Q2m *= 0.5;
90
91    maxQ2[kk] = std::min(limitQ2, Q2m);
92    TableCrossSec[ONQ2*kk] = 0.0;
93
94//    G4cout<<" kk  eGeV[kk] "<<kk<<"  "<<eGeV[kk]<<G4endl;
95  }
96}
97
98/////////////////////////////////////////////////////////////////////////
99//
100//
101
102void G4ElasticData::DefineNucleusParameters(G4double A)
103{
104  switch (AtomicWeight)
105    {
106    case 207:
107    case 208:
108      R1       = 20.5;
109//      R1       = 17.5;
110//      R1       = 21.3;   
111      R2       = 15.74;
112//      R2       = 10.74;
113
114      Pnucl    = 0.4;
115      Aeff     = 0.7;
116      break;
117    case 237:
118    case 238:
119      R1       = 21.7;   
120      R2       = 16.5;
121      Pnucl    = 0.4;
122      Aeff     = 0.7;
123      break;
124    case 90:
125    case 91:
126      R1    = 16.5*1.0;
127      R2    = 11.62;
128      Pnucl = 0.4;
129      Aeff  = 0.7;
130      break;
131    case 58:
132    case 59:
133      R1    = 15.0*1.05;
134      R2    = 9.9;
135      Pnucl = 0.45;
136      Aeff  = 0.85;
137      break;
138    case 48:
139    case 47:
140      R1    = 14.0;
141      R2    = 9.26;
142      Pnucl = 0.31;
143      Aeff  = 0.75;
144      break;
145    case 40:
146    case 41:
147      R1    = 13.3;
148      R2    = 9.26;
149      Pnucl = 0.31;
150      Aeff  = 0.75;
151      break;
152    case 28:
153    case 29:
154      R1    = 12.0;
155      R2    = 7.64;
156      Pnucl = 0.253;
157      Aeff  = 0.8;
158      break;
159    case 16:
160      R1    = 10.50;
161      R2    = 5.5;
162      Pnucl = 0.7;
163      Aeff  = 0.98;
164      break;
165    case 12:
166      R1    = 9.3936;
167      R2    = 4.63;
168      Pnucl = 0.7;
169//      Pnucl = 0.5397;
170      Aeff  = 1.0;
171      break;
172    case 11:
173      R1    = 9.0;
174      R2    = 5.42;
175      Pnucl = 0.19;
176//      Pnucl = 0.5397;
177      Aeff  = 0.9;
178      break;
179    case 9:
180      R1    = 9.9;
181      R2    = 6.5;
182      Pnucl = 0.690;
183      Aeff  = 0.95;
184      break;
185    case 4:
186      R1    = 5.3;   
187      R2    = 3.7;
188      Pnucl = 0.4;
189      Aeff  = 0.75;
190      break;
191    case 1:
192      R1    = 4.5;   
193      R2    = 2.3;
194      Pnucl = 0.177;
195      Aeff  = 0.9;
196      break;
197    default:
198      R1    = 4.45*std::pow(A - 1.,0.309)*0.9;
199      R2    = 2.3 *std::pow(A, 0.36);
200
201      if(A < 100 && A > 3) Pnucl = 0.176 + 0.00275*A;
202      else                 Pnucl = 0.4;
203
204//G4cout<<" Deault: A= "<<A<<"  R1 R2 Aeff Pnucl "<<R1<<"  "<<R2<<"  "
205//      <<Aeff<<"  "<<Pnucl<<G4endl;
206
207      if(A >= 100)               Aeff = 0.7;
208      else if(A < 100 && A > 75) Aeff = 1.5 - 0.008*A;
209      else                       Aeff = 0.9;
210      break;
211    }
212//G4cout<<" Result: A= "<<A<<"  R1 R2 Aeff Pnucl "<<R1<<"  "<<R2<<"  "
213//      <<Aeff<<"  "<<Pnucl<<G4endl;
214}
215
216////////////////////////////////////////////////////////////////////
217//
218//  The constructor for the generating of events
219//
220
221G4ElasticHadrNucleusHE::G4ElasticHadrNucleusHE()
222  :G4HadronicInteraction("G4ElasticHadrNucleusHE")
223{
224  verboseLevel = 0;
225  plabLowLimit = 20.0*MeV;
226  lowestEnergyLimit = 0.0;
227
228  MbToGeV2  =  2.568;
229  sqMbToGeV =  1.602;
230  Fm2ToGeV2 =  25.68;
231  GeV2      =  GeV*GeV;
232  protonM   =  proton_mass_c2/GeV;
233  protonM2  =  protonM*protonM;
234
235   BoundaryP[0]=9.0;BoundaryTG[0]=5.0;BoundaryTL[0]=0.;
236   BoundaryP[1]=20.0;BoundaryTG[1]=1.5;BoundaryTL[1]=0.;
237   BoundaryP[2]=5.0; BoundaryTG[2]=1.0;BoundaryTL[2]=1.5;
238   BoundaryP[3]=8.0; BoundaryTG[3]=3.0;BoundaryTL[3]=0.;
239   BoundaryP[4]=7.0; BoundaryTG[4]=3.0;BoundaryTL[4]=0.;
240   BoundaryP[5]=5.0; BoundaryTG[5]=2.0;BoundaryTL[5]=0.;
241   BoundaryP[6]=5.0; BoundaryTG[6]=1.5;BoundaryTL[6]=3.0;
242
243  Binom();
244  // energy in GeV
245  Energy[0] = 0.4;
246  Energy[1] = 0.6;
247  Energy[2] = 0.8;
248  LowEdgeEnergy[0] = 0.0;
249  LowEdgeEnergy[1] = 0.5;
250  LowEdgeEnergy[2] = 0.7;
251  G4double e = 1.0;
252  G4double f = std::pow(10.,0.1);
253  for(G4int i=3; i<NENERGY; i++) {
254    Energy[i] = e;
255    LowEdgeEnergy[i] = e/f;
256    e *= f*f;
257  }
258  nistManager = G4NistManager::Instance();
259
260  // PDG code for hadrons
261  G4int ic[NHADRONS] = {211,-211,2112,2212,321,-321,130,310,311,-311,
262                        3122,3222,3112,3212,3312,3322,3334,
263                        -2212,-2112,-3122,-3222,-3112,-3212,-3312,-3322,-3334};
264  // internal index
265  G4int id[NHADRONS] = {2,3,6,0,4,5,4,4,4,5,
266                        0,0,0,0,0,0,0,
267                        1,7,1,1,1,1,1,1,1};
268
269  G4int id1[NHADRONS] = {3,4,1,0,5,6,5,5,5,6,
270                        0,0,0,0,0,0,0,
271                        2,2,2,2,2,2,2,2,2};
272
273  for(G4int j=0; j<NHADRONS; j++) 
274  {
275    HadronCode[j]  = ic[j];
276    HadronType[j]  = id[j];
277    HadronType1[j] = id1[j];
278
279    for(G4int k = 0; k < 93; k++) SetOfElasticData[j][k] = 0;
280  } 
281}
282
283///////////////////////////////////////////////////////////////////
284//
285//
286
287G4ElasticHadrNucleusHE::~G4ElasticHadrNucleusHE()
288{
289  for(G4int j = 0; j < NHADRONS; j++) 
290  {
291    for(G4int k = 0; k < 93; k++) 
292    {
293      if(SetOfElasticData[j][k]) delete SetOfElasticData[j][k];
294    }
295  } 
296}
297
298////////////////////////////////////////////////////////////////////
299//
300//
301
302G4HadFinalState * G4ElasticHadrNucleusHE::ApplyYourself(
303                          const  G4HadProjectile  &aTrack,
304                                 G4Nucleus        &targetNucleus)
305{
306  theParticleChange.Clear();
307
308  const G4HadProjectile* aParticle = &aTrack;
309  G4double ekin = aParticle->GetKineticEnergy();
310
311  if( ekin <= lowestEnergyLimit ) 
312  {
313    theParticleChange.SetEnergyChange(ekin);
314    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
315    return &theParticleChange;
316  }
317
318  G4double aTarget = targetNucleus.GetN();
319  G4double zTarget = targetNucleus.GetZ();
320
321  G4double plab = aParticle->GetTotalMomentum();
322
323  if (verboseLevel >1)
324  { 
325    G4cout << "G4ElasticHadrNucleusHE: Incident particle plab=" 
326           << plab/GeV << " GeV/c " 
327           << " ekin(MeV) = " << ekin/MeV << "  " 
328           << aParticle->GetDefinition()->GetParticleName() << G4endl;
329  }
330  // Scattered particle referred to axis of incident particle
331
332  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
333  G4double m1 = theParticle->GetPDGMass();
334
335  G4int Z = static_cast<G4int>(zTarget+0.5);
336  G4int A = static_cast<G4int>(aTarget+0.5);
337  G4int projPDG = theParticle->GetPDGEncoding();
338
339  if (verboseLevel>1) 
340  {
341    G4cout << "G4ElasticHadrNucleusHE for " << theParticle->GetParticleName()
342           << " PDGcode= " << projPDG << " on nucleus Z= " << Z
343           << " A= " << A
344           << G4endl;
345  }
346  G4ParticleDefinition * theDef = 0;
347
348  if      (Z == 1 && A == 1) theDef = G4Proton::Proton();
349  else if (Z == 1 && A == 2) theDef = G4Deuteron::Deuteron();
350  else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
351  else if (Z == 2 && A == 3) theDef = G4He3::He3();
352  else if (Z == 2 && A == 4) theDef = G4Alpha::Alpha();
353  else                       theDef = G4ParticleTable::
354                               GetParticleTable()->FindIon(Z,A,0,Z);
355 
356  G4double m2 = theDef->GetPDGMass();
357  G4LorentzVector lv1 = aParticle->Get4Momentum();
358  G4LorentzVector lv(0.0,0.0,0.0,m2);   
359  lv += lv1;
360
361  G4ThreeVector bst = lv.boostVector();
362  lv1.boost(-bst);
363
364  G4ThreeVector p1 = lv1.vect();
365  G4double ptot = p1.mag();
366  G4double tmax = 4.0*ptot*ptot;
367  G4double t = 0.0;
368
369  // Choose generator
370  G4bool swave = false;
371
372  // S-wave for very low energy
373  if(plab < plabLowLimit) swave = true;
374
375  // normal sampling
376  if(!swave) {
377    t = SampleT(theParticle,plab,Z,A);
378    if(t > tmax) swave = true;
379  }
380
381  if(swave) t = G4UniformRand()*tmax;
382
383  // Sampling in CM system
384  G4double phi  = G4UniformRand()*twopi;
385  G4double cost = 1. - 2.0*t/tmax;
386  G4double sint;
387
388  if( cost >= 1.0 ) 
389  {
390    cost = 1.0;
391    sint = 0.0;
392  }
393  else if( cost <= -1.0) 
394  {
395    cost = -1.0;
396    sint =  0.0;
397  }
398  else 
399  {
400    sint = std::sqrt((1.0-cost)*(1.0+cost));
401  } 
402  if (verboseLevel>1)
403    G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
404 
405  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
406  v1 *= ptot;
407  G4LorentzVector nlv1( v1.x(), v1.y(), v1.z(), std::sqrt(ptot*ptot + m1*m1));
408
409  nlv1.boost(bst); 
410
411  G4double eFinal = nlv1.e() - m1;
412
413  if (verboseLevel > 1)
414  { 
415    G4cout << "Scattered: "
416           << nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal
417           << " Proj: 4-mom " << lv1
418           <<G4endl;
419  }
420  if( eFinal < 0.0 ) 
421  {
422    G4cout << "G4ElasticHadrNucleusHE WARNING ekin= " << eFinal
423           << " after scattering of " 
424           << aParticle->GetDefinition()->GetParticleName()
425           << " p(GeV/c)= " << plab
426           << " on " << theDef->GetParticleName()
427           << G4endl;
428    eFinal = 0.0;
429    nlv1.setE(m1);
430  }
431
432  theParticleChange.SetMomentumChange( nlv1.vect().unit() );
433  theParticleChange.SetEnergyChange(eFinal);
434 
435  G4LorentzVector nlv0 = lv - nlv1;
436  G4double erec =  nlv0.e() - m2;
437
438  if (verboseLevel > 1)
439  { 
440    G4cout << "Recoil: "
441           << nlv0<<" m= " << m2 << " ekin(MeV)= " << erec
442           <<G4endl;
443  }
444  if(erec < 0.0) 
445  {
446    G4cout << "G4ElasticHadrNucleusHE WARNING Erecoil(MeV)= " << erec
447           << " after scattering of " 
448           << aParticle->GetDefinition()->GetParticleName()
449           << " p(GeV/c)= " << plab
450           << " on " << theDef->GetParticleName()
451           << G4endl;
452    nlv0.setE(m2);
453  }
454  G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
455  theParticleChange.AddSecondary(aSec);
456
457  return &theParticleChange;
458}
459
460//////////////////////////////////////////////////////////////////////////
461//
462//
463
464G4double G4ElasticHadrNucleusHE::
465                     SampleT(  const G4ParticleDefinition* p,
466                                     G4double inLabMom, G4int Z, G4int N)
467{
468  G4double plab  = inLabMom/GeV;   // (GeV/c)
469  G4double Q2 = 0;
470
471  iHadrCode = p->GetPDGEncoding();
472
473  NumbN = N;
474
475  if(verboseLevel > 1)
476  {
477    G4cout<< " G4ElasticHadrNucleusHE::SampleT: " 
478          << " for " << p->GetParticleName() 
479          << " at Z= " << Z << " A= " << N
480          << " plab(GeV)= " << plab
481          << G4endl;
482  }
483  G4int idx;
484
485  for( idx = 0 ; idx < NHADRONS; idx++) 
486  {
487    if(iHadrCode == HadronCode[idx]) break;
488  }
489
490  // Hadron is not in the list
491
492  if( idx >= NHADRONS ) return Q2;
493
494  iHadron = HadronType[idx];
495  iHadrCode = HadronCode[idx];
496
497  if(Z==1)
498    {
499      hMass  = p->GetPDGMass()/GeV;
500      hMass2 = hMass*hMass;
501
502      G4double T = sqrt(plab*plab+hMass2)-hMass;
503
504      if(T > 0.4) Q2 = HadronProtonQ2(p, plab);
505
506      if (verboseLevel>1)
507        G4cout<<"  Proton : Q2  "<<Q2<<G4endl;
508    }
509  else
510    {
511      G4ElasticData* ElD1 = SetOfElasticData[idx][Z];
512
513      // Construct elastic data
514      if(!ElD1) 
515        {
516          G4double AWeight = nistManager->GetAtomicMassAmu(Z);
517          ElD1 = new  G4ElasticData(p, Z, AWeight, Energy);
518          SetOfElasticData[idx][Z] = ElD1;
519   
520          if(verboseLevel > 1)
521            {
522              G4cout<< " G4ElasticHadrNucleusHE::SampleT:  new record " << idx
523                    << " for " << p->GetParticleName() << " Z= " << Z
524                    << G4endl;
525            }
526        } 
527      hMass          = ElD1->massGeV;
528      hMass2         = ElD1->mass2GeV2;
529      G4double M     = ElD1->massA;
530      G4double M2    = ElD1->massA2;
531      G4double plab2 = plab*plab;
532      G4double Q2max = 4.*plab2*M2/
533        (hMass2 + M2 + 2.*M*std::sqrt(plab2 + hMass2));
534
535      // sample scattering
536      G4double T = sqrt(plab2+hMass2)-hMass;
537
538      if(T > 0.4) Q2 = HadronNucleusQ2_2(ElD1, Z, plab, Q2max);
539
540      if(verboseLevel > 1)
541        G4cout<<" SampleT: Q2(GeV^2)= "<<Q2<< "  t/tmax= " << Q2/Q2max <<G4endl;
542    }
543  return  Q2*GeV2;
544}
545
546//////////////////////////////////////////////////////////////////////////
547//
548//
549
550G4double G4ElasticHadrNucleusHE::
551                          HadronNucleusQ2_2(G4ElasticData* pElD, G4int Z, 
552                                            G4double plab, G4double tmax)
553{
554  G4double LineFq2[ONQ2];
555
556  G4double Rand = G4UniformRand();
557
558  G4int      iNumbQ2 = 0;
559  G4double   Q2 = 0.0;
560
561  G4double ptot2 = plab*plab;
562  G4double ekin  = std::sqrt(hMass2 + ptot2) - hMass;
563
564  if(verboseLevel > 1)
565    G4cout<<"Q2_2: ekin  plab  "<<ekin<<"    "<<plab<<"  tmax "<<tmax<<G4endl;
566
567  // Find closest energy bin
568  G4int NumbOnE; 
569  for( NumbOnE = 0; NumbOnE < NENERGY-1; NumbOnE++ ) 
570  {
571    if( ekin <= LowEdgeEnergy[NumbOnE+1] ) break;
572  }
573  G4double* dNumbQ2 = pElD->TableQ2;
574
575  G4int index = NumbOnE*ONQ2;
576
577  // Select kinematics for node energy
578  G4double T     = Energy[NumbOnE];
579  hLabMomentum2  = T*(T + 2.*hMass);
580  G4double Q2max = pElD->maxQ2[NumbOnE];
581  G4int length   = pElD->dnkE[NumbOnE];
582
583  // Build vector
584  if(length == 0) 
585    {
586      R1    = pElD->R1;
587      R2    = pElD->R2;
588      Aeff  = pElD->Aeff;
589      Pnucl = pElD->Pnucl;
590      hLabMomentum = std::sqrt(hLabMomentum2);
591 
592      DefineHadronValues(Z);
593
594      if(verboseLevel >0)
595        {
596          G4cout<<"1  plab  T  "<<plab<<"  "<<T<<"  sigTot  B  ReIm  "
597                <<HadrTot<<"  "<<HadrSlope<<"  "<<HadrReIm<<G4endl;
598          G4cout<<"  R1  R2  Aeff  p  "<<R1<<"  "<<R2<<"  "<<Aeff<<"  "
599                <<Pnucl<<G4endl;
600        }
601
602      pElD->CrossSecMaxQ2[NumbOnE] = 1.0;
603
604      if(verboseLevel > 1)
605        G4cout<<" HadrNucleusQ2_2: NumbOnE= " << NumbOnE
606              << " length= " << length
607              << " Q2max= " << Q2max
608              << " ekin= " << ekin <<G4endl;
609   
610      pElD->TableCrossSec[index] = 0;
611
612
613      dQ2 = pElD->TableQ2[1]-pElD->TableQ2[0];
614
615      GetHeavyFq2(NumbN, LineFq2);  //  %%%%%%%%%%%%%%%%%%%%%%%%%
616
617      for(G4int ii=0; ii<ONQ2; ii++)
618        {
619          //if(verboseLevel > 2)
620          //  G4cout<<"  ii LineFq2  "<<ii<<"  "<<LineFq2[ii]/LineFq2[ONQ2-1]
621          //    <<"  dF(q2) "<<LineFq2[ii]-LineFq2[ii-1]<<G4endl;
622
623          pElD->TableCrossSec[index+ii] = LineFq2[ii]/LineFq2[ONQ2-1];
624        }
625   
626      pElD->dnkE[NumbOnE] = ONQ2;
627      length = ONQ2;
628    } 
629
630  G4double* dNumbFQ2 = &(pElD->TableCrossSec[index]);
631
632  for( iNumbQ2 = 1; iNumbQ2<length; iNumbQ2++ ) 
633    {
634      if(Rand <= pElD->TableCrossSec[index+iNumbQ2]) break;
635    }
636  Q2 = GetQ2_2(iNumbQ2, dNumbQ2, dNumbFQ2, Rand);
637
638  if(tmax < Q2max) Q2 *= tmax/Q2max;
639
640  if(verboseLevel > 1)
641    G4cout<<" HadrNucleusQ2_2(2): Q2= "<<Q2<<" iNumbQ2= " << iNumbQ2
642          << " rand= " << Rand << G4endl;
643 
644  return Q2;
645}       
646
647///////////////////////////////////////////////////////////////////////
648//
649//  The randomization of one dimensional array
650//
651G4double G4ElasticHadrNucleusHE::GetQ2_2(G4int kk, G4double * Q,
652                                         G4double * F, G4double ranUni)
653{
654  G4double ranQ2;
655
656  G4double F1  = F[kk-2];
657  G4double F2  = F[kk-1];
658  G4double F3  = F[kk];
659  G4double X1  = Q[kk-2];
660  G4double X2  = Q[kk-1];
661  G4double X3  = Q[kk];
662
663  if(verboseLevel > 2) 
664    G4cout << "GetQ2_2 kk= " << kk << " X2= " << X2 << " X3= " << X3
665           << " F2= " << F2 << " F3= " << F3 << " Rndm= " << ranUni << G4endl;
666
667  if(kk == 1 || kk == 0)
668  {
669     F1  = F[0]; 
670     F2  = F[1];
671     F3  = F[2];
672     X1  = Q[0];
673     X2  = Q[1];
674     X3  = Q[2];
675  }
676
677  G4double F12 = F1*F1;
678  G4double F22 = F2*F2;
679  G4double F32 = F3*F3;
680
681  G4double D0  = F12*F2+F1*F32+F3*F22-F32*F2-F22*F1-F12*F3;
682
683  if(verboseLevel > 2) 
684    G4cout << "       X1= " << X1 << " F1= " << F1 << "  D0= " 
685           << D0 << G4endl; 
686
687  if(std::abs(D0) < 0.00000001)
688    { 
689      ranQ2 = X2 + (ranUni - F2)*(X3 - X2)/(F3 - F2);
690    }
691  else   
692    {
693      G4double DA = X1*F2+X3*F1+X2*F3-X3*F2-X1*F3-X2*F1;
694      G4double DB = X2*F12+X1*F32+X3*F22-X2*F32-X3*F12-X1*F22;
695      G4double DC = X3*F2*F12+X2*F1*F32+X1*F3*F22
696                   -X1*F2*F32-X2*F3*F12-X3*F1*F22;
697      ranQ2 = (DA*ranUni*ranUni + DB*ranUni + DC)/D0;
698    }
699  return ranQ2;         //  MeV^2
700}
701
702////////////////////////////////////////////////////////////////////////
703//
704//
705G4double G4ElasticHadrNucleusHE::GetHeavyFq2(G4int Nucleus, G4double* LineF)
706{
707  G4int ii, jj, aSimp;
708  G4double curQ2, curSec;
709  G4double curSum = 0.0;
710  G4double totSum = 0.0;
711
712  G4double ddQ2 = dQ2/20;
713  G4double Q2l  = 0;
714
715  LineF[0] = 0;
716  for(ii = 1; ii<ONQ2; ii++)
717    {
718      curSum = 0;
719      aSimp  = 4;   
720
721      for(jj = 0; jj<20; jj++)
722        {
723          curQ2 = Q2l+jj*ddQ2;
724
725          curSec  = HadrNucDifferCrSec(Nucleus, curQ2);
726          curSum += curSec*aSimp;
727
728          if(aSimp > 3) aSimp = 2;
729          else          aSimp = 4;
730
731          if(jj == 0 && verboseLevel>2)
732            G4cout<<"  Q2  "<<curQ2<<"  AIm  "<<aAIm<<"  DIm  "<<aDIm
733                  <<"  Diff  "<<curSec<<"  totSum  "<<totSum<<G4endl;
734        }
735
736      Q2l    += dQ2;
737      curSum *= ddQ2/2.3;   //  $$$$$$$$$$$$$$$$$$$$$$$
738      totSum += curSum;
739
740      LineF[ii] = totSum;
741       
742      if (verboseLevel>2)
743        G4cout<<"  GetHeavy: Q2  dQ2  totSum  "<<Q2l<<"  "<<dQ2<<"  "<<totSum
744              <<"  curSec  "
745              <<curSec<<"  totSum  "<< totSum<<"  DTot "
746              <<curSum<<G4endl;
747    }     
748  return totSum;
749}
750
751////////////////////////////////////////////////////////////////////////
752//
753//
754
755G4double G4ElasticHadrNucleusHE::GetLightFq2(G4int Z, G4int Nucleus, 
756                                             G4double Q2)
757{
758  // Scattering of proton
759  if(Z == 1) 
760  {
761    G4double SqrQ2  = std::sqrt(Q2);
762    G4double ConstU = 2.*(hMass2 + protonM2) - Q2;
763
764    G4double y = (1.-Coeff1-Coeff0)/HadrSlope*(1.-std::exp(-HadrSlope*Q2))
765      + Coeff0*(1.-std::exp(-Slope0*Q2))
766      + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1.)
767      + 2.*Coeff1/Slope1*(1./Slope1-(1./Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
768
769    return y;
770  }
771
772  // The preparing of probability function 
773
774  G4double prec = Nucleus > 208  ?  1.0e-7 : 1.0e-6;
775
776  G4double    Stot     = HadrTot*MbToGeV2;     //  Gev^-2
777  G4double    Bhad     = HadrSlope;         //  GeV^-2
778  G4double    Asq      = 1+HadrReIm*HadrReIm;
779  G4double    Rho2     = std::sqrt(Asq);
780
781//  if(verboseLevel >1)
782    G4cout<<" Fq2 Before for i Tot B Im "<<HadrTot<<"  "<<HadrSlope<<"  "
783      <<HadrReIm<<G4endl;
784
785  if(verboseLevel > 1) {
786    G4cout << "GetFq2: Stot= " << Stot << " Bhad= " << Bhad
787           <<"  Im "<<HadrReIm
788           << " Asq= " << Asq << G4endl;
789    G4cout << "R1= " << R1 << " R2= " << R2 << " Pnucl= " << Pnucl <<G4endl;
790  }
791  G4double    R12      = R1*R1;
792  G4double    R22      = R2*R2;
793  G4double    R12B     = R12+2*Bhad;
794  G4double    R22B     = R22+2*Bhad;
795
796  G4double    Norm     = (R12*R1-Pnucl*R22*R2); // HP->Aeff;
797
798  G4double    R13      = R12*R1/R12B;
799  G4double    R23      = Pnucl*R22*R2/R22B;
800  G4double    Unucl    = Stot/twopi/Norm*R13;
801  G4double    UnucRho2 = -Unucl*Rho2;
802
803  G4double    FiH      = std::asin(HadrReIm/Rho2);
804  G4double    NN2      = R23/R13;
805
806  if(verboseLevel > 2) 
807    G4cout << "UnucRho2= " << UnucRho2 << " FiH= " << FiH << " NN2= " << NN2
808           << " Norm= " << Norm << G4endl;
809
810  G4double    dddd;
811 
812  G4double    Prod0    = 0;
813  G4double    N1       = -1.0;
814  G4double    Tot0     = 0;
815  G4double    exp1;
816
817  G4double    Prod3 ;
818  G4double    exp2  ;
819  G4double    N4, N5, N2, Prod1, Prod2;
820  G4int    i1, i2, m1, m2;
821
822  for(i1 = 1; i1<= Nucleus; i1++) ////++++++++++  i1
823    {
824      N1    = -N1*Unucl*(Nucleus-i1+1)/i1*Rho2;
825      Prod1 = 0;
826      Tot0  = 0;
827      N2    = -1;
828
829      for(i2 = 1; i2<=Nucleus; i2++) ////+++++++++ i2
830        {
831          N2    = -N2*Unucl*(Nucleus-i2+1)/i2*Rho2;
832          Prod2 = 0; 
833          N5    = -1/NN2;
834          for(m2=0; m2<= i2; m2++) ////+++++++++ m2
835            {
836              Prod3 = 0;
837              exp2  = 1/(m2/R22B+(i2-m2)/R12B);
838              N5    = -N5*NN2;
839              N4    = -1/NN2;
840              for(m1=0; m1<=i1; m1++) ////++++++++ m1
841                {
842                  exp1  = 1/(m1/R22B+(i1-m1)/R12B);
843                  dddd  = exp1+exp2;
844                  N4    = -N4*NN2;
845                  Prod3 = Prod3+N4*exp1*exp2*
846                    (1-std::exp(-Q2*dddd/4))/dddd*4*SetBinom[i1][m1];
847               }                                   // m1
848              Prod2 = Prod2 +Prod3*N5*SetBinom[i2][m2];
849            }                                      // m2
850          Prod1 = Prod1 + Prod2*N2*std::cos(FiH*(i1-i2));
851
852          if (std::fabs(Prod2*N2/Prod1)<prec) break;
853        }                                         // i2
854      Prod0   = Prod0 + Prod1*N1;
855      if(std::fabs(N1*Prod1/Prod0) < prec) break;
856
857    }                                           // i1
858
859  Prod0 *= 0.25*pi/MbToGeV2;  //  This is in mb
860
861  if(verboseLevel>1) 
862    G4cout << "GetLightFq2 Z= " << Z << " A= " << Nucleus
863           <<" Q2= " << Q2 << " Res= " << Prod0 << G4endl;
864  return Prod0;
865}
866//  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
867G4double G4ElasticHadrNucleusHE::
868   HadrNucDifferCrSec(G4int Nucleus, G4double aQ2)
869{
870//   ------ All external kinematical variables are in MeV -------
871//            ------ but internal in GeV !!!!  ------
872
873  G4double    theQ2 = aQ2;   ///GeV/GeV; 
874
875  // Scattering of proton
876  if(Nucleus == 1) 
877  {
878    G4double SqrQ2  = std::sqrt(aQ2);
879    G4double ConstU = hMass2 + protonM2-2*protonM*HadrEnergy - aQ2;
880
881    G4double MaxT = 4*MomentumCM*MomentumCM;
882
883     BoundaryTL[0] = MaxT;
884     BoundaryTL[1] = MaxT;
885     BoundaryTL[3] = MaxT;
886     BoundaryTL[4] = MaxT;
887     BoundaryTL[5] = MaxT;
888
889    G4double dSigPodT;
890
891    dSigPodT = HadrTot*HadrTot*(1+HadrReIm*HadrReIm)*
892                 (
893                  Coeff1*std::exp(-Slope1*SqrQ2)+
894                  Coeff2*std::exp( Slope2*(ConstU)+aQ2)+
895                  (1-Coeff1-Coeff0)*std::exp(-HadrSlope*aQ2)+
896                 +Coeff0*std::exp(-Slope0*aQ2)
897//                +0.1*(1-std::fabs(CosTh))
898                  )/16/3.1416*2.568;
899
900    return dSigPodT;
901  }
902
903    G4double    Stot     = HadrTot*MbToGeV2; 
904    G4double    Bhad     = HadrSlope; 
905    G4double    Asq      = 1+HadrReIm*HadrReIm;
906    G4double    Rho2     = std::sqrt(Asq);
907    G4double    Pnuclp   = 0.001;
908                Pnuclp   = Pnucl;
909    G4double    R12      = R1*R1;
910    G4double    R22      = R2*R2;
911    G4double    R12B     = R12+2*Bhad;
912    G4double    R22B     = R22+2*Bhad;
913    G4double    R12Ap    = R12+20;
914    G4double    R22Ap    = R22+20;
915    G4double    R13Ap    = R12*R1/R12Ap;
916    G4double    R23Ap    = R22*R2/R22Ap*Pnuclp;
917    G4double    R23dR13  = R23Ap/R13Ap;
918    G4double    R12Apd   = 2/R12Ap;
919    G4double    R22Apd   = 2/R22Ap;
920    G4double R12ApdR22Ap = 0.5*(R12Apd+R22Apd);
921
922    G4double DDSec1p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R1/4));
923    G4double DDSec2p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/
924                             std::sqrt((R12+R22)/2)/4));
925    G4double DDSec3p  = (DDSect2+DDSect3*std::log(1.06*2*HadrEnergy/R2/4));
926
927    G4double    Norm     = (R12*R1-Pnucl*R22*R2)*Aeff;
928    G4double    Normp    = (R12*R1-Pnuclp*R22*R2)*Aeff;
929    G4double    R13      = R12*R1/R12B;
930    G4double    R23      = Pnucl*R22*R2/R22B;
931    G4double    Unucl    = Stot/twopi/Norm*R13;
932    G4double    UnuclScr = Stot/twopi/Normp*R13Ap;
933    G4double    SinFi    = HadrReIm/Rho2;
934    G4double    FiH      = std::asin(SinFi);
935    G4double    N        = -1;
936    G4double    N2       = R23/R13;
937
938    G4double    ImElasticAmpl0  = 0;
939    G4double    ReElasticAmpl0  = 0;
940
941    G4double    exp1;
942    G4double    N4;
943    G4double    Prod1, Tot1=0, medTot, DTot1, DmedTot;
944    G4int       i;
945
946    for( i=1; i<=Nucleus; i++)
947    {
948      N       = -N*Unucl*(Nucleus-i+1)/i*Rho2;
949      N4      = 1;
950      Prod1   = std::exp(-theQ2/i*R12B/4)/i*R12B;
951      medTot  = R12B/i;
952
953       for(G4int l=1; l<=i; l++)
954       {
955         exp1    = l/R22B+(i-l)/R12B;
956         N4      = -N4*(i-l+1)/l*N2;
957         Prod1   = Prod1+N4/exp1*std::exp(-theQ2/exp1/4);
958         medTot  = medTot+N4/exp1;
959       }  // end l
960
961      ReElasticAmpl0  = ReElasticAmpl0+Prod1*N*std::sin(FiH*i);
962      ImElasticAmpl0  = ImElasticAmpl0+Prod1*N*std::cos(FiH*i);
963      Tot1            = Tot1+medTot*N*std::cos(FiH*i);
964      if(std::fabs(Prod1*N/ImElasticAmpl0) < 0.000001) break;
965    }      // i
966
967    ImElasticAmpl0 = ImElasticAmpl0*pi/2.568;   // The amplitude in mB
968    ReElasticAmpl0 = ReElasticAmpl0*pi/2.568;   // The amplitude in mB
969    Tot1           = Tot1*twopi/2.568;
970
971    G4double C1 = R13Ap*R13Ap/2*DDSec1p;
972    G4double C2 = 2*R23Ap*R13Ap/2*DDSec2p;
973    G4double C3 = R23Ap*R23Ap/2*DDSec3p;
974
975    G4double N1p  = 1;
976
977    G4double Din1 = 0.5;     
978
979    Din1  = 0.5*(C1*std::exp(-theQ2/8*R12Ap)/2*R12Ap-
980                 C2/R12ApdR22Ap*std::exp(-theQ2/4/R12ApdR22Ap)+
981                 C3*R22Ap/2*std::exp(-theQ2/8*R22Ap));
982
983    DTot1 = 0.5*(C1/2*R12Ap-C2/R12ApdR22Ap+C3*R22Ap/2);
984
985    G4double exp1p;
986    G4double exp2p;
987    G4double exp3p;
988    G4double N2p;
989    G4double Din2, BinCoeff;
990
991    BinCoeff = 1;
992
993    for( i = 1; i<= Nucleus-2; i++)
994    {
995      N1p     = -N1p*UnuclScr*(Nucleus-i-1)/i*Rho2;
996      N2p     = 1;
997      Din2    = 0;
998      DmedTot = 0;
999        for(G4int l = 0; l<=i; l++) 
1000        {
1001          if(l == 0)      BinCoeff = 1;
1002          else if(l !=0 ) BinCoeff = BinCoeff*(i-l+1)/l;
1003
1004          exp1  = l/R22B+(i-l)/R12B;
1005          exp1p = exp1+R12Apd;
1006          exp2p = exp1+R12ApdR22Ap;
1007          exp3p = exp1+R22Apd;
1008
1009          Din2  = Din2 + N2p*BinCoeff*
1010                  (C1/exp1p*std::exp(-theQ2/4/exp1p)-
1011                   C2/exp2p*std::exp(-theQ2/4/exp2p)+
1012                   C3/exp3p*std::exp(-theQ2/4/exp3p));
1013
1014          DmedTot = DmedTot + N2p*BinCoeff*
1015            (C1/exp1p-C2/exp2p+C3/exp3p);
1016
1017          N2p   = -N2p*R23dR13;
1018        }     // l
1019
1020        Din1  = Din1+Din2*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
1021        DTot1 = DTot1+DmedTot*N1p/*Mnoj[i]*//(i+2)/(i+1)*std::cos(FiH*i);
1022 
1023        if(std::fabs(Din2*N1p/Din1) < 0.000001) break;
1024    }           //  i
1025
1026    Din1 = -Din1*Nucleus*(Nucleus-1)
1027                 /2/pi/Normp/2/pi/Normp*16*pi*pi;
1028
1029    DTot1 = DTot1*Nucleus*(Nucleus-1)
1030                 /2/pi/Normp/2/pi/Normp*16*pi*pi;
1031
1032    DTot1 *= 5;   //  $$$$$$$$$$$$$$$$$$$$$$$$
1033//     Din1  *= 0.2;    //   %%%%%%%%%%%%%%%%%%%%%%%   proton
1034//     Din1 *= 0.05;    //   %%%%%%%%%%%%%%%%%%%%%%%  pi+
1035//  ----------------  dSigma/d|-t|,  mb/(GeV/c)^-2  -----------------
1036
1037    G4double DiffCrSec2 = (ReElasticAmpl0*ReElasticAmpl0+
1038                           (ImElasticAmpl0+Din1)*
1039                           (ImElasticAmpl0+Din1))*2/4/pi;
1040
1041    Tot1   = Tot1-DTot1;
1042     //  Tott1  = Tot1*1.0;
1043    Dtot11 = DTot1;
1044    aAIm   = ImElasticAmpl0;
1045    aDIm   = Din1;
1046
1047    return DiffCrSec2*1.0;  //  dSig/d|-t|,  mb/(GeV/c)^-2
1048}   // function
1049//  ##############################################
1050
1051////////////////////////////////////////////////////////////////
1052//
1053//
1054
1055void  G4ElasticHadrNucleusHE::DefineHadronValues(G4int Z)
1056{
1057  HadrEnergy = std::sqrt(hMass2 + hLabMomentum2);
1058
1059  G4double sHadr = 2.*HadrEnergy*protonM+protonM2+hMass2;
1060  G4double sqrS  = std::sqrt(sHadr);
1061  G4double Ecm   = 0.5*(sHadr-hMass2+protonM2)/sqrS;
1062  MomentumCM     = std::sqrt(Ecm*Ecm-protonM2);
1063 
1064  if(verboseLevel>2)
1065    G4cout << "GetHadrVall.: Z= " << Z << " iHadr= " << iHadron
1066           << " E(GeV)= " << HadrEnergy << " sqrS= " << sqrS
1067           << " plab= " << hLabMomentum   
1068           <<"  E - m  "<<HadrEnergy - hMass<< G4endl;
1069
1070  G4double TotN = 0.0;
1071  G4double logE = std::log(HadrEnergy);
1072  G4double logS = std::log(sHadr);
1073           TotP = 0.0;
1074
1075  switch (iHadron)
1076    {
1077    case 0:                  //  proton, neutron
1078    case 6:
1079
1080      if(hLabMomentum > 10)
1081        TotP = TotN = 7.5*logE - 40.12525 + 103*std::pow(sHadr,-0.165); //  mb
1082
1083      else
1084        {
1085// ==================  neutron  ================
1086
1087////      if(iHadrCode == 2112)
1088
1089
1090          if( hLabMomentum > 1.4 )
1091            TotN = 33.3+15.2*(hLabMomentum2-1.35)/
1092              (std::pow(hLabMomentum,2.37)+0.95);
1093               
1094          else if(hLabMomentum > 0.8)
1095            {
1096              G4double A0 = logE + 0.0513;
1097              TotN = 33.0 + 25.5*A0*A0; 
1098            }
1099          else 
1100            {
1101              G4double A0 = logE - 0.2634;  // log(1.3)
1102              TotN = 33.0 + 30.*A0*A0*A0*A0;
1103            }
1104//  =================  proton  ===============
1105//       else if(iHadrCode == 2212)
1106          {
1107            if(hLabMomentum >= 1.05)
1108              {
1109                TotP = 39.0+75.*(hLabMomentum-1.2)/
1110                  (hLabMomentum2*hLabMomentum+0.15);
1111              }
1112
1113            else if(hLabMomentum >= 0.7)
1114              {
1115                 G4double A0 = logE + 0.3147;
1116                 TotP = 23.0 + 40.*A0*A0;
1117              }
1118            else 
1119              {
1120                TotP = 23.+50.*std::pow(std::log(0.73/hLabMomentum),3.5);
1121              }
1122          }
1123        }
1124
1125//      HadrTot = 0.5*(82*TotP+126*TotN)/104;  //  $$$$$$$$$$$$$$$$$$
1126      HadrTot = 0.5*(TotP+TotN);
1127//  ...................................................
1128      //  Proton slope
1129      if(hLabMomentum >= 2.)       HadrSlope = 5.44 + 0.88*logS;
1130
1131      else if(hLabMomentum >= 0.5) HadrSlope = 3.73*hLabMomentum-0.37;
1132
1133      else                         HadrSlope = 1.5;
1134
1135//  ...................................................
1136      if(hLabMomentum >= 1.2)
1137         HadrReIm  = 0.13*(logS - 5.8579332)*std::pow(sHadr,-0.18);
1138       
1139      else if(hLabMomentum >= 0.6) 
1140        HadrReIm = -75.5*(std::pow(hLabMomentum,0.25)-0.95)/
1141          (std::pow(3*hLabMomentum,2.2)+1);     
1142
1143      else 
1144        HadrReIm = 15.5*hLabMomentum/(27*hLabMomentum2*hLabMomentum+2);
1145//  ...................................................
1146      DDSect2   = 2.2;                              //mb*GeV-2
1147      DDSect3   = 0.6;                               //mb*GeV-2
1148      //  ================== lambda  ==================
1149      if( iHadrCode == 3122)
1150        {
1151          HadrTot   *= 0.88;
1152          HadrSlope *=0.85;
1153        }
1154      //  ================== sigma +  ==================
1155      else if( iHadrCode == 3222)
1156        {
1157          HadrTot   *=0.81;
1158          HadrSlope *=0.85;
1159        }
1160      //  ================== sigma 0,-  ==================
1161      else if(iHadrCode == 3112 || iHadrCode == 3212 )
1162        {
1163          HadrTot   *=0.88;
1164          HadrSlope *=0.85;
1165        }
1166      //  ===================  xi  =================
1167      else if( iHadrCode == 3312 || iHadrCode == 3322 )
1168        {
1169          HadrTot   *=0.77;
1170          HadrSlope *=0.75;
1171        }
1172      //  =================  omega  =================
1173      else if( iHadrCode == 3334)
1174        {
1175          HadrTot   *=0.78;
1176          HadrSlope *=0.7;
1177        }
1178
1179      break;
1180//  ===========================================================
1181    case 1:              //   antiproton
1182    case 7:              //   antineutron
1183
1184      HadrTot   = 5.2+5.2*logE + 123.2/sqrS;     //  mb
1185      HadrSlope = 8.32+0.57*logS;                //(GeV/c)^-2
1186
1187      if( HadrEnergy < 1000 )
1188        HadrReIm  = 0.06*(sqrS-2.236)*(sqrS-14.14)*std::pow(sHadr,-0.8);
1189      else
1190        HadrReIm  = 0.6*(logS - 5.8579332)*std::pow(sHadr,-0.25);
1191
1192      DDSect2   = 11;                            //mb*(GeV/c)^-2
1193      DDSect3   = 3;                             //mb*(GeV/c)^-2
1194      //  ================== lambda  ==================
1195      if( iHadrCode == -3122)
1196        {
1197          HadrTot   *= 0.88;
1198          HadrSlope *=0.85;
1199        }
1200      //  ================== sigma +  ==================
1201      else if( iHadrCode == -3222)
1202        {
1203          HadrTot   *=0.81;
1204          HadrSlope *=0.85;
1205        }
1206      //  ================== sigma 0,-  ==================
1207      else if(iHadrCode == -3112 || iHadrCode == -3212 )
1208        {
1209          HadrTot   *=0.88;
1210          HadrSlope *=0.85;
1211        }
1212      //  ===================  xi  =================
1213      else if( iHadrCode == -3312 || iHadrCode == -3322 )
1214        {
1215          HadrTot   *=0.77;
1216          HadrSlope *=0.75;
1217        }
1218      //  =================  omega  =================
1219      else if( iHadrCode == -3334)
1220        {
1221          HadrTot   *=0.78;
1222          HadrSlope *=0.7;
1223        }
1224
1225      break;
1226//  -------------------------------------------
1227    case 2:             //   pi plus, pi minus
1228    case 3:
1229
1230      if(hLabMomentum >= 3.5)
1231        TotP = 10.6+2.*logE + 25.*std::pow(HadrEnergy,-0.43); // mb
1232//  =========================================
1233      else if(hLabMomentum >= 1.15)
1234        {
1235          G4double x = (hLabMomentum - 2.55)/0.55; 
1236          G4double y = (hLabMomentum - 1.47)/0.225;
1237          TotP = 3.2*std::exp(-x*x) + 12.*std::exp(-y*y) + 27.5;
1238        }
1239//  =========================================
1240      else if(hLabMomentum >= 0.4)
1241        {
1242        TotP  = 88*(logE+0.2877)*(logE+0.2877)+14.0;
1243        }
1244//  =========================================
1245      else 
1246        {
1247          G4double x = (hLabMomentum - 0.29)/0.085;
1248          TotP = 20. + 180.*std::exp(-x*x);
1249        }
1250//  -------------------------------------------
1251
1252      if(hLabMomentum >= 3.0 )
1253        TotN = 10.6 + 2.*logE + 30.*std::pow(HadrEnergy,-0.43); // mb
1254
1255      else if(hLabMomentum >= 1.3) 
1256        {
1257          G4double x = (hLabMomentum - 2.1)/0.4;
1258          G4double y = (hLabMomentum - 1.4)/0.12;
1259          TotN = 36.1+0.079 - 4.313*logE + 3.*std::exp(-x*x) + 
1260                                              1.5*std::exp(-y*y);
1261        }
1262      else if(hLabMomentum >= 0.65)
1263        {
1264          G4double x = (hLabMomentum - 0.72)/0.06;
1265          G4double y = (hLabMomentum - 1.015)/0.075;
1266          TotN = 36.1 + 10.*std::exp(-x*x) + 24*std::exp(-y*y);
1267        }
1268      else if(hLabMomentum >= 0.37)
1269        {
1270          G4double x = std::log(hLabMomentum/0.48);
1271          TotN = 26. + 110.*x*x;
1272        }
1273      else 
1274        {
1275          G4double x = (hLabMomentum - 0.29)/0.07;
1276          TotN = 28.0 + 40.*std::exp(-x*x);
1277        }
1278      HadrTot = (TotP+TotN)/2;
1279//  ........................................
1280      HadrSlope = 7.28+0.245*logS;        // GeV-2
1281      HadrReIm  = 0.2*(logS - 4.6051702)*std::pow(sHadr,-0.15);
1282
1283      DDSect2   = 0.7;                               //mb*GeV-2
1284      DDSect3   = 0.27;                              //mb*GeV-2
1285
1286      break;
1287//  ==========================================================
1288    case 4:            //  K plus
1289
1290      HadrTot   = 10.6+1.8*logE + 9.0*std::pow(HadrEnergy,-0.55);  // mb
1291      if(HadrEnergy>100) HadrSlope = 15.0;
1292      else HadrSlope = 1.0+1.76*logS - 2.84/sqrS;   // GeV-2
1293
1294      HadrReIm  = 0.4*(sHadr-20)*(sHadr-150)*std::pow(sHadr+50,-2.1);
1295      DDSect2   = 0.7;                             //mb*GeV-2
1296      DDSect3   = 0.21;                            //mb*GeV-2
1297      break;
1298//  =========================================================
1299     case 5:              //   K minus
1300
1301       HadrTot   = 10+1.8*logE + 25./sqrS; // mb
1302       HadrSlope = 6.98+0.127*logS;         // GeV-2
1303       HadrReIm  = 0.4*(sHadr-20)*(sHadr-20)*std::pow(sHadr+50,-2.1);
1304       DDSect2   = 0.7;                             //mb*GeV-2
1305       DDSect3   = 0.27;                            //mb*GeV-2
1306       break;
1307  }   
1308//  =========================================================
1309  if(verboseLevel>2)
1310    G4cout << "HadrTot= " << HadrTot << " HadrSlope= " << HadrSlope
1311           << " HadrReIm= " << HadrReIm << " DDSect2= " << DDSect2
1312           << " DDSect3= " << DDSect3 << G4endl;
1313
1314  if(Z != 1) return;
1315
1316  // Scattering of protons
1317
1318  Coeff0 = Coeff1 = Coeff2 = 0.0;
1319  Slope0 = Slope1 = 1.0;
1320  Slope2 = 5.0;
1321
1322  // data for iHadron=0
1323  static const G4double EnP0[6]={1.5,3.0,5.0,9.0,14.0,19.0};
1324  static const G4double C0P0[6]={0.15,0.02,0.06,0.08,0.0003,0.0002};
1325  static const G4double C1P0[6]={0.05,0.02,0.03,0.025,0.0,0.0};
1326  static const G4double B0P0[6]={1.5,2.5,3.0,4.5,1.4,1.25};
1327  static const G4double B1P0[6]={5.0,1.0,3.5,4.0,4.8,4.8};
1328     
1329  // data for iHadron=6,7
1330  static const G4double EnN[5]={1.5,5.0,10.0,14.0,20.0};
1331  static const G4double C0N[5]={0.0,0.0,0.02,0.02,0.01};
1332  static const G4double C1N[5]={0.06,0.008,0.0015,0.001,0.0003};
1333  static const G4double B0N[5]={1.5,2.5,3.8,3.8,3.5};
1334  static const G4double B1N[5]={1.5,2.2,3.6,4.5,4.8};
1335
1336  // data for iHadron=1
1337  static const G4double EnP[2]={1.5,4.0};
1338  static const G4double C0P[2]={0.001,0.0005};
1339  static const G4double C1P[2]={0.003,0.001};
1340  static const G4double B0P[2]={2.5,4.5};
1341  static const G4double B1P[2]={1.0,4.0};
1342
1343  // data for iHadron=2
1344  static const G4double EnPP[4]={1.0,2.0,3.0,4.0};
1345  static const G4double C0PP[4]={0.0,0.0,0.0,0.0};
1346  static const G4double C1PP[4]={0.15,0.08,0.02,0.01};
1347  static const G4double B0PP[4]={1.5,2.8,3.8,3.8};
1348  static const G4double B1PP[4]={0.8,1.6,3.6,4.6};
1349
1350  // data for iHadron=3
1351  static const G4double EnPPN[4]={1.0,2.0,3.0,4.0};
1352  static const G4double C0PPN[4]={0.0,0.0,0.0,0.0};
1353  static const G4double C1PPN[4]={0.0,0.0,0.0,0.0};
1354  static const G4double B0PPN[4]={1.5,2.8,3.8,3.8};
1355  static const G4double B1PPN[4]={0.8,1.6,3.6,4.6};
1356
1357  // data for iHadron=4
1358  static const G4double EnK[4]={1.4,2.33,3.0,5.0};
1359  static const G4double C0K[4]={0.0,0.0,0.0,0.0};
1360  static const G4double C1K[4]={0.01,0.007,0.005,0.003};
1361  static const G4double B0K[4]={1.5,2.0,3.8,3.8};
1362  static const G4double B1K[4]={1.6,1.6,1.6,1.6};
1363
1364  // data for iHadron=5
1365  static const G4double EnKM[2]={1.4,4.0};
1366  static const G4double C0KM[2]={0.006,0.002};
1367  static const G4double C1KM[2]={0.00,0.00};
1368  static const G4double B0KM[2]={2.5,3.5};
1369  static const G4double B1KM[2]={1.6,1.6};
1370
1371  switch(iHadron)
1372    {
1373    case 0 :
1374
1375      if(hLabMomentum <BoundaryP[0])
1376        InterpolateHN(6,EnP0,C0P0,C1P0,B0P0,B1P0);
1377
1378      Coeff2 = 0.8/hLabMomentum/hLabMomentum;
1379      break; 
1380
1381    case  6 :
1382
1383      if(hLabMomentum < BoundaryP[1])
1384        InterpolateHN(5,EnN,C0N,C1N,B0N,B1N);
1385
1386      Coeff2 = 0.8/hLabMomentum/hLabMomentum;
1387      break; 
1388
1389    case 1 :
1390    case 7 :
1391      if(hLabMomentum <  BoundaryP[2])
1392        InterpolateHN(2,EnP,C0P,C1P,B0P,B1P);
1393      break; 
1394
1395    case 2 :
1396
1397      if(hLabMomentum < BoundaryP[3])
1398        InterpolateHN(4,EnPP,C0PP,C1PP,B0PP,B1PP);
1399
1400      Coeff2 = 0.02/hLabMomentum;
1401      break; 
1402
1403    case 3 :
1404
1405      if(hLabMomentum < BoundaryP[4])
1406        InterpolateHN(4,EnPPN,C0PPN,C1PPN,B0PPN,B1PPN);
1407
1408      Coeff2 = 0.02/hLabMomentum;
1409      break;
1410 
1411    case 4 :
1412
1413      if(hLabMomentum < BoundaryP[5])
1414        InterpolateHN(4,EnK,C0K,C1K,B0K,B1K);
1415
1416      if(hLabMomentum < 1) Coeff2 = 0.34;
1417      else  Coeff2 = 0.34/hLabMomentum2/hLabMomentum;
1418      break; 
1419
1420    case 5 :
1421      if(hLabMomentum < BoundaryP[6])
1422        InterpolateHN(2,EnKM,C0KM,C1KM,B0KM,B1KM);
1423
1424      if(hLabMomentum < 1) Coeff2 = 0.01;
1425      else  Coeff2 = 0.01/hLabMomentum2/hLabMomentum;
1426      break; 
1427    }
1428
1429  if(verboseLevel > 2) 
1430    G4cout<<"  HadrVal : Plasb  "<<hLabMomentum
1431          <<"  iHadron  "<<iHadron<<"  HadrTot  "<<HadrTot<<G4endl;
1432}
1433
1434//  =====================================================
1435void  G4ElasticHadrNucleusHE::
1436       GetKinematics(const G4ParticleDefinition * aHadron,
1437                           G4double MomentumH)
1438{
1439  if (verboseLevel>1)
1440    G4cout<<"1  GetKin.: HadronName MomentumH "
1441          <<aHadron->GetParticleName()<<"  "<<MomentumH<<G4endl;
1442
1443  DefineHadronValues(1);
1444
1445  G4double Sh     = 2.0*protonM*HadrEnergy+protonM2+hMass2; // GeV
1446
1447  ConstU = 2*protonM2+2*hMass2-Sh;
1448
1449  G4double MaxT = 4*MomentumCM*MomentumCM;
1450
1451  BoundaryTL[0] = MaxT; //2.0;
1452  BoundaryTL[1] = MaxT;
1453  BoundaryTL[3] = MaxT;
1454  BoundaryTL[4] = MaxT;
1455  BoundaryTL[5] = MaxT;
1456
1457  G4int NumberH=0;
1458
1459  while(iHadrCode!=HadronCode[NumberH]) NumberH++;
1460
1461  NumberH = HadronType1[NumberH];   
1462
1463  if(MomentumH<BoundaryP[NumberH]) MaxTR = BoundaryTL[NumberH];
1464  else MaxTR = BoundaryTG[NumberH];
1465
1466  if (verboseLevel>1)
1467    G4cout<<"3  GetKin. : NumberH  "<<NumberH
1468          <<"  Bound.P[NumberH] "<<BoundaryP[NumberH]
1469          <<"  Bound.TL[NumberH] "<<BoundaryTL[NumberH]
1470          <<"  Bound.TG[NumberH] "<<BoundaryTG[NumberH]
1471          <<"  MaxT MaxTR "<<MaxT<<"  "<<MaxTR<<G4endl;
1472
1473//     GetParametersHP(aHadron, MomentumH);
1474}
1475//  ============================================================
1476G4double G4ElasticHadrNucleusHE::GetFt(G4double Q2)
1477{
1478  G4double Fdistr=0;
1479  G4double SqrQ2 = std::sqrt(Q2);
1480 
1481  Fdistr = (1-Coeff1-Coeff0) //-0.0*Coeff2*std::exp(ConstU))
1482    /HadrSlope*(1-std::exp(-HadrSlope*Q2))
1483    + Coeff0*(1-std::exp(-Slope0*Q2))
1484    + Coeff2/Slope2*std::exp(Slope2*ConstU)*(std::exp(Slope2*Q2)-1)
1485    + 2*Coeff1/Slope1*(1/Slope1-(1/Slope1+SqrQ2)*std::exp(-Slope1*SqrQ2));
1486
1487  if (verboseLevel>1)
1488    G4cout<<"Old:  Coeff0 Coeff1 Coeff2 "<<Coeff0<<"  "
1489          <<Coeff1<<"  "<<Coeff2<<"  Slope Slope0 Slope1 Slope2 "
1490          <<HadrSlope<<"  "<<Slope0<<"  "<<Slope1<<"  "<<Slope2
1491          <<"  Fdistr "<<Fdistr<<G4endl;
1492  return Fdistr;
1493}
1494//  +++++++++++++++++++++++++++++++++++++++
1495G4double G4ElasticHadrNucleusHE::GetQ2(G4double Ran)
1496{
1497  G4double DDD0=MaxTR*0.5, DDD1=0.0, DDD2=MaxTR, delta;
1498  G4double Q2=0;
1499
1500  FmaxT = GetFt(MaxTR);
1501  delta = GetDistrFun(DDD0)-Ran;
1502
1503  while(std::fabs(delta) > 0.0001)
1504    {
1505      if(delta>0) 
1506      {
1507        DDD2 = DDD0;
1508        DDD0 = (DDD0+DDD1)*0.5;
1509      }
1510      else if(delta<0)
1511      {
1512        DDD1 = DDD0; 
1513        DDD0 = (DDD0+DDD2)*0.5;
1514      }
1515      delta = GetDistrFun(DDD0)-Ran;
1516    }
1517
1518  Q2 = DDD0;
1519
1520  return Q2;
1521}
1522//  ++++++++++++++++++++++++++++++++++++++++++
1523G4double G4ElasticHadrNucleusHE::
1524       HadronProtonQ2(const G4ParticleDefinition * p,
1525                            G4double inLabMom)
1526{
1527
1528  hMass         = p->GetPDGMass()/GeV;
1529  hMass2        = hMass*hMass;
1530  hLabMomentum  = inLabMom;
1531  hLabMomentum2 = hLabMomentum*hLabMomentum;
1532  HadrEnergy    = sqrt(hLabMomentum2+hMass2); 
1533
1534  G4double Rand = G4UniformRand();
1535
1536  GetKinematics(p, inLabMom);
1537
1538  G4double Q2 = GetQ2(Rand);
1539
1540  return Q2;
1541}
1542
1543//  ===========================================
1544
1545///////////////////////////////////////////////////////////////////
1546//
1547// 
1548
1549void G4ElasticHadrNucleusHE::Binom()
1550{
1551  G4int N, M;
1552  G4double  Fact1, J;
1553
1554  for(N = 0; N < 240; N++)
1555  {
1556    J = 1;
1557
1558      for( M = 0; M <= N; M++ )
1559      {
1560          Fact1 = 1;
1561
1562          if ( ( N > 0 ) && ( N > M ) && M > 0 )
1563          {
1564              J *= G4double(N-M+1)/G4double(M);
1565              Fact1 = J;
1566          }
1567          SetBinom[N][M] = Fact1;
1568      }
1569  }
1570  return;
1571}
1572
1573
1574//
1575//
1576///////////////////////////////////////////////////////////
1577
Note: See TracBrowser for help on using the repository browser.