ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLEtaNElasticChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLEtaNElasticChannel.cc
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 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 #include "G4INCLEtaNElasticChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 
45 namespace G4INCL {
46 
47  EtaNElasticChannel::EtaNElasticChannel(Particle *p1, Particle *p2)
48  : particle1(p1), particle2(p2)
49  {
50 
51  }
52 
53  EtaNElasticChannel::~EtaNElasticChannel(){
54 
55  }
56 
57  void EtaNElasticChannel::fillFinalState(FinalState *fs) {
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  }
208 }