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

Last change on this file since 1358 was 1340, checked in by garnier, 14 years ago

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