Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4INCL::EtaNElasticChannel Class Reference

#include <G4INCLEtaNElasticChannel.hh>

Inheritance diagram for G4INCL::EtaNElasticChannel:
Inheritance graph
[legend]
Collaboration diagram for G4INCL::EtaNElasticChannel:
Collaboration graph
[legend]

Public Member Functions

 EtaNElasticChannel (Particle *, Particle *)
 
virtual ~EtaNElasticChannel ()
 
void fillFinalState (FinalState *fs)
 
- Public Member Functions inherited from G4INCL::IChannel
 IChannel ()
 
virtual ~IChannel ()
 
FinalState * getFinalState ()
 

Detailed Description

Definition at line 47 of file G4INCLEtaNElasticChannel.hh.

Constructor & Destructor Documentation

G4INCL::EtaNElasticChannel::EtaNElasticChannel ( Particle *  p1,
Particle *  p2 
)

Definition at line 47 of file G4INCLEtaNElasticChannel.cc.

48  : particle1(p1), particle2(p2)
49  {
50 
51  }
G4INCL::EtaNElasticChannel::~EtaNElasticChannel ( )
virtual

Definition at line 53 of file G4INCLEtaNElasticChannel.cc.

53  {
54 
55  }

Member Function Documentation

void G4INCL::EtaNElasticChannel::fillFinalState ( FinalState *  fs)
virtual

Implements G4INCL::IChannel.

Definition at line 57 of file G4INCLEtaNElasticChannel.cc.

57  {
58  Particle * nucleon;
59  Particle * eta;
60  if(particle1->isNucleon()) {
61  nucleon = particle1;
62  eta = particle2;
63  } else {
64  nucleon = particle2;
65  eta = particle1;
66  }
67 
68  G4double plab=KinematicsUtils::momentumInLab(particle1, particle2);
69 
70  G4double sh=nucleon->getEnergy()+eta->getEnergy();
71  G4double mn=nucleon->getMass();
72  G4double me=eta->getMass();
73  G4double en=(sh*sh+mn*mn-me*me)/(2*sh);
74  nucleon->setEnergy(en);
75  G4double ee=std::sqrt(en*en-mn*mn+me*me);
76  eta->setEnergy(ee);
77  G4double pn=std::sqrt(en*en-mn*mn);
78 
79  ThreeVector mom_nucleon;
80 
81  if (plab < 250.) {
82 // Isotropy
83  mom_nucleon = Random::normVector(pn);
84  }
85 
86 // From Kamano
87  else {
88 
89  const G4double pi=std::acos(-1.0);
90  G4double x1;
91  G4double u1;
92  G4double fteta;
93  G4double teta;
94  G4double fi;
95 
96  G4double a0;
97  G4double a1;
98  G4double a2;
99  G4double a3;
100  G4double a4;
101  G4double a5;
102  G4double a6;
103 
104  if (plab > 1400.) plab=1400.; // no information on angular distributions above plab=1400 MeV
105  G4double p6=std::pow(plab, 6);
106  G4double p5=std::pow(plab, 5);
107  G4double p4=std::pow(plab, 4);
108  G4double p3=std::pow(plab, 3);
109  G4double p2=std::pow(plab, 2);
110  G4double p1=plab;
111 
112 // a6
113  if (plab < 300.) {
114  a6=-8.384000E-08*p1 - 1.15452E-04;
115  }
116  else if (plab < 500.){
117  a6=1.593966E-13*p4 - 2.619560E-10*p3 + 1.564701E-07*p2 - 3.986627E-05*p1 + 3.622575E-03;
118  }
119  else {
120  a6=6.143615E-20*p6 - 3.157181E-16*p5 + 6.348289E-13*p4 - 6.117961E-10*p3 + 2.764542E-07*p2 - 4.391048E-05*p1 - 1.443857E-03;
121  }
122 // a5
123  if (plab < 650.) {
124  a5=-9.021076E-18*p6 + 2.176771E-14*p5 - 2.136095E-11*p4 + 1.100580E-08*p3 - 3.150857E-06*p2 + 4.761016E-04*p1 - 2.969608E-02;
125  }
126  else if (plab < 950.){
127  a5=4.424756E-18*p6 - 1.756295E-14*p5 + 2.625428E-11*p4 - 1.678272E-08*p3 + 2.227237E-06*p2 + 2.146666E-03*p1 - 7.065712E-01;
128  }
129  else {
130  a5=2.209585E-19*p6 - 1.546647E-15*p5 + 4.578142E-12*p4 - 7.303856E-09*p3 + 6.604074E-06*p2 - 3.205628E-03*p1 + 6.534893E-01;
131  }
132 // a4
133  if (plab < 700.) {
134  a4=4.826684E-17*p6 - 1.534471E-13*p5 + 1.907868E-10*p4 - 1.192317E-07*p3 + 3.988902E-05*p2 - 6.822100E-03*p1 + 4.684685E-01;
135  }
136  else {
137  a4=-3.245143E-18*p6 + 2.174395E-14*p5 - 6.012288E-11*p4 + 8.772790E-08*p3 - 7.113554E-05*p2 + 3.029285E-02*p1 - 5.237677E+00;
138  }
139 // a3
140  if (plab < 650.) {
141  a3=3.783071E-17*p6 - 1.151454E-13*p5 + 1.357165E-10*p4 - 8.036891E-08*p3 + 2.572396E-05*p2 - 4.245566E-03*p1 + 2.832772E-01;
142  }
143  else {
144  a3=-5.063316E-18*p6 + 3.223757E-14*p5 - 8.435635E-11*p4 + 1.159487E-07*p3 - 8.812510E-05*p2 + 3.500692E-02*p1 - 5.624556E+00;
145  }
146 // a2
147  if (plab < 500.) {
148  a2=-6.085067E-14*p5 + 1.354078E-10*p4 - 1.124158E-07*p3 + 4.292106E-05*p2 - 7.218145E-03*p1 + 4.584962E-01;
149  }
150  else if (plab < 750.) {
151  a2= 9.512730E-11*p4 - 2.362724E-07*p3 + 2.171883E-04*p2 - 8.742722E-02*p1 + 1.309433E+01;
152  }
153  else {
154  a2=-4.228889E-18*p6 + 2.798222E-14*p5 - 7.640831E-11*p4 + 1.100124E-07*p3 - 8.778573E-05*p2 + 3.652772E-02*p1 - 6.025497E+00;
155  }
156 // a1
157  if (plab < 500.) {
158  a1=-1.524408E-14*p5 + 3.007021E-11*p4 - 2.129570E-08*p3 + 5.607250E-06*p2 - 3.001598E-04*p1 + 8.701280E-04;
159  }
160  else if (plab < 750.) {
161  a1=-3.255396E-11*p4 + 8.168681E-08*p3 - 7.447474E-05*p2 + 2.917630E-02*p1 - 4.152037E+00;
162  }
163  else {
164  a1=9.964504E-19*p6 - 6.380168E-15*p5 + 1.638691E-11*p4 - 2.107063E-08*p3 + 1.347462E-05*p2 - 3.318304E-03*p1 - 5.030932E-02;
165  }
166 // a0
167  a0=-3.220143E-17*p6 + 1.789654E-13*p5 - 3.912863E-10*p4 + 4.181510E-07*p3 - 2.147259E-04*p2 + 3.856266E-02*p1 + 2.609971E+00;
168 
169  G4double interg1=2.*(a6/7. + a4/5. + a2/3. + a0); // (integral to normalize)
170  G4double f1=(a6+a5+a4+a3+a2+a1+a0)/interg1; // (Max normalized)
171 
172  G4int passe1=0;
173  while (passe1==0) {
174  // Sample x from -1 to 1
175  x1=Random::shoot();
176  if (Random::shoot() > 0.5) x1=-x1;
177 
178  // Sample u from 0 to 1
179  u1=Random::shoot();
180  fteta=(a6*x1*x1*x1*x1*x1*x1+a5*x1*x1*x1*x1*x1+a4*x1*x1*x1*x1+a3*x1*x1*x1+a2*x1*x1+a1*x1+a0)/interg1;
181  // The condition
182  if (u1*f1 < fteta) {
183  teta=std::acos(x1);
184  // std::cout << x1 << " " << fteta << " "<< f1/interg1 << " " << u1 << " " << interg1 << std::endl;
185  passe1=1;
186  }
187  }
188 
189  fi=(2.0*pi)*Random::shoot();
190 
191  ThreeVector mom_nucleon1(
192  pn*std::sin(teta)*std::cos(fi),
193  pn*std::sin(teta)*std::sin(fi),
194  pn*std::cos(teta)
195  );
196 
197  mom_nucleon = -mom_nucleon1 ;
198 
199  }
200 
201  nucleon->setMomentum(mom_nucleon);
202  eta->setMomentum(-mom_nucleon);
203 
204  fs->addModifiedParticle(nucleon);
205  fs->addModifiedParticle(eta);
206 
207  }
a0
const G4double a0
Definition: G4IonCoulombCrossSection.cc:68
G4int
int G4int
Definition: G4Types.hh:78
G4INCL::Random::normVector
ThreeVector normVector(G4double norm=1.)
Definition: G4INCLRandom.cc:134
G4InuclParticleNames::nucleon
G4bool nucleon(G4int ityp)
Definition: G4InuclParticleNames.hh:95
G4INCL::Particle::isNucleon
G4bool isNucleon() const
Definition: G4INCLParticle.hh:232
G4INCL::Random::shoot
G4double shoot()
Definition: G4INCLRandom.cc:93
pi
static constexpr double pi
Definition: G4SIunits.hh:75
G4INCL::KinematicsUtils::momentumInLab
G4double momentumInLab(Particle const *const p1, Particle const *const p2)
gives the momentum in the lab frame of two particles.
Definition: G4INCLKinematicsUtils.cc:134
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

The documentation for this class was generated from the following files: