ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLEtaNToPiNChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLEtaNToPiNChannel.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 "G4INCLEtaNToPiNChannel.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  EtaNToPiNChannel::EtaNToPiNChannel(Particle *p1, Particle *p2)
48  : particle1(p1), particle2(p2)
49  {
50 
51  }
52 
53  EtaNToPiNChannel::~EtaNToPiNChannel(){
54 
55  }
56 
57  void EtaNToPiNChannel::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  const G4double r2 = Random::shoot();
71  if (nucleon->getType() == Neutron) {
72  if (r2*3. < 2.) {
73  nucleon->setType(Proton);
74  eta->setType(PiMinus);
75  }
76  else {
77  nucleon->setType(Neutron);
78  eta->setType(PiZero);
79  }
80  }
81  else {
82  if (r2*3. < 2.) {
83  nucleon->setType(Neutron);
84  eta->setType(PiPlus);
85  }
86  else {
87  nucleon->setType(Proton);
88  eta->setType(PiZero);
89  }
90  }
91 
92  G4double sh=nucleon->getEnergy()+eta->getEnergy();
93  G4double mn=nucleon->getMass();
94  G4double me=eta->getMass();
95  G4double en=(sh*sh+mn*mn-me*me)/(2*sh);
96  nucleon->setEnergy(en);
97  G4double ee=std::sqrt(en*en-mn*mn+me*me);
98  eta->setEnergy(ee);
99  G4double pn=std::sqrt(en*en-mn*mn);
100 
101  const G4double pi=std::acos(-1.0);
102  G4double x1;
103  G4double u1;
104  G4double fteta;
105  G4double teta;
106  G4double fi;
107 
108  G4double a0;
109  G4double a1;
110  G4double a2;
111  G4double a3;
112  G4double a4;
113  G4double a5;
114  G4double a6;
115 
116  if (plab > 1400.) plab=1400.; // no information on angular distributions above plab=1400 MeV
117  G4double p6=std::pow(plab, 6);
118  G4double p5=std::pow(plab, 5);
119  G4double p4=std::pow(plab, 4);
120  G4double p3=std::pow(plab, 3);
121  G4double p2=std::pow(plab, 2);
122  G4double p1=plab;
123 
124  // a6
125  if (plab <= 600.) {
126  a6=5.721872E-18*p6 - 1.063594E-14*p5 +
127  7.812226E-12*p4 - 2.947343E-09*p3 +
128  5.955500E-07*p2 - 6.081534E-05*p1 + 2.418893E-03;
129  }
130  else {
131  a6=1.549323E-18*p6 - 9.570613E-15*p5 +
132  2.428560E-11*p4 - 3.237490E-08*p3 +
133  2.385312E-05*p2 - 9.167580E-03*p1 + 1.426952E+00;
134  }
135  // a5
136  if (plab <= 700.) {
137  a5=-3.858406E-16*p6 + 7.397533E-13*p5 -
138  5.344420E-10*p4 + 1.865842E-07*p3 -
139  3.234292E-05*p2 + 2.552380E-03*p1 - 6.810842E-02;
140  }
141  else {
142  a5=-3.775268E-17*p6 + 2.445059E-13*p5 -
143  6.503137E-10*p4 + 9.065678E-07*p3 -
144  6.953576E-04*p2 + 2.757524E-01*p1 - 4.328028E+01;
145  }
146  // a4
147  if (plab <= 550.) {
148  a4=-2.051840E-16*p6 + 3.858551E-13*p5 -
149  3.166229E-10*p4 + 1.353545E-07*p3 -
150  2.631251E-05*p2 + 2.109593E-03*p1 - 5.633076E-02;
151  }
152  else if (plab <= 650.) {
153  a4=-1.698136E-05*p2 + 1.827203E-02*p1 - 4.482122E+00;
154  }
155  else {
156  a4=-2.808337E-17*p6 + 1.640033E-13*p5 -
157  3.820460E-10*p4 + 4.452787E-07*p3 -
158  2.621981E-04*p2 + 6.530743E-02*p1 - 2.447717E+00;
159  }
160  // a3
161  if (plab <= 700.) {
162  a3=7.061866E-16*p6 - 1.356389E-12*p5 +
163  9.783322E-10*p4 - 3.407333E-07*p3 +
164  5.903545E-05*p2 - 4.735559E-03*p1 + 1.270435E-01;
165  }
166  else {
167  a3=1.138088E-16*p6 - 7.459580E-13*p5 +
168  2.015156E-09*p4 - 2.867416E-06*p3 +
169  2.261028E-03*p2 - 9.323442E-01*p1 + 1.552846E+02;
170  }
171  // a2
172  if (plab <= 550.) {
173  a2=1.352952E-17*p6 - 3.030435E-13*p5 +
174  4.624668E-10*p4 - 2.759605E-07*p3 +
175  6.996373E-05*p2 - 4.745692E-03*p1 + 1.524349E-01;
176  }
177  else if (plab <= 700.) {
178  a2=5.514651E-08*p3 - 8.734112E-05*p2 + 4.108704E-02*p1 - 5.116601E+00;
179  }
180  else {
181  a2=5.621795E-17*p6 - 3.701960E-13*p5 +
182  1.005796E-09*p4 - 1.441294E-06*p3 +
183  1.146234E-03*p2 - 4.775194E-01*p1 + 8.084776E+01;
184  }
185  // a1
186  if (plab <= 500.) {
187  a1=-2.425827E-16*p6 + 4.113350E-13*p5 -
188  2.342298E-10*p4 + 4.934322E-08*p3 -
189  3.564530E-06*p2 + 6.516398E-04*p1 + 2.547230E-01;
190  }
191  else if (plab <= 700.) {
192  a1=-1.824213E-10*p4 + 3.599251E-07*p3 -
193  2.480862E-04*p2 + 6.894931E-02*p1 - 5.760562E+00;
194  }
195  else {
196  a1=-5.139366E-17*p6 + 3.408224E-13*p5 -
197  9.341903E-10*p4 + 1.354028E-06*p3 -
198  1.093509E-03*p2 + 4.653326E-01*p1 - 8.068436E+01;
199  }
200  // a0
201  if (plab <= 400.) {
202  a0=1.160837E-13*p6 - 1.813002E-10*p5 +
203  1.155391E-07*p4 - 3.862737E-05*p3 +
204  7.230513E-03*p2 - 7.469799E-01*p1 + 3.830064E+01;
205  }
206  else if (plab <= 700.) {
207  a0=2.267918E-14*p6 - 7.593899E-11*p5 +
208  1.049849E-07*p4 - 7.669301E-05*p3 +
209  3.123846E-02*p2 - 6.737221E+00*p1 + 6.032010E+02;
210  }
211  else {
212  a0=-1.851188E-17*p6 + 1.281122E-13*p5 -
213  3.686161E-10*p4 + 5.644116E-07*p3 -
214  4.845757E-04*p2 + 2.203918E-01*p1 - 4.100383E+01;
215  }
216 
217  G4double interg1=2.*(a6/7. + a4/5. + a2/3. + a0); // (integral to normalize)
218  G4double f1=(a6+a5+a4+a3+a2+a1+a0)/interg1; // (Max normalized)
219 
220  G4int passe1=0;
221  while (passe1==0) {
222  // Sample x from -1 to 1
223  x1=Random::shoot();
224  if (Random::shoot() > 0.5) x1=-x1;
225 
226  // Sample u from 0 to 1
227  u1=Random::shoot();
228  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;
229  // The condition
230  if (u1*f1 < fteta) {
231  teta=std::acos(x1);
232  // std::cout << x1 << " " << fteta << " "<< f1/interg1 << " " << u1 << " " << interg1 << std::endl;
233  passe1=1;
234  }
235  }
236 
237  fi=(2.0*pi)*Random::shoot();
238 
239  ThreeVector mom_nucleon(
240  pn*std::sin(teta)*std::cos(fi),
241  pn*std::sin(teta)*std::sin(fi),
242  pn*std::cos(teta)
243  );
244  // end real distribution
245 
246  nucleon->setMomentum(-mom_nucleon);
247  eta->setMomentum(mom_nucleon);
248 
249  fs->addModifiedParticle(nucleon);
250  fs->addModifiedParticle(eta);
251  }
252 
253 }