ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLElasticChannel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLElasticChannel.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 "G4INCLElasticChannel.hh"
39 #include "G4INCLRandom.hh"
40 #include "G4INCLKinematicsUtils.hh"
41 #include "G4INCLParticleTable.hh"
42 #include "G4INCLCrossSections.hh"
43 #include "G4INCLGlobals.hh"
44 
45 namespace G4INCL {
46 
48  :particle1(p1), particle2(p2)
49  {
50  }
51 
53  {
54  }
55 
57  {
58  ParticleType p1TypeOld = particle1->getType();
59  ParticleType p2TypeOld = particle2->getType();
60 
61  /* Concerning the way we calculate the lab momentum, see the considerations
62  * in CrossSections::elasticNNLegacy().
63  */
66 
67  const G4int isospin = ParticleTable::getIsospin(particle1->getType()) +
69 
70  // Calculate the outcome of the channel:
72  G4double pnorm = std::sqrt(psq);
74  G4double btmax = 4.0 * psq * b;
75  G4double z = std::exp(-btmax);
76  G4double ranres = Random::shoot();
77  G4double y = 1.0 - ranres * (1.0 - z);
78  G4double T = std::log(y)/b;
79  G4int iexpi = 0;
80  G4double apt = 1.0;
81 
82  // Handle np case
83  if((particle1->getType() == Proton && particle2->getType() == Neutron) ||
84  (particle1->getType() == Neutron && particle2->getType() == Proton)) {
85  if(pl > 800.0) {
86  const G4double x = 0.001 * pl; // Transform to GeV
87  apt = (800.0/pl)*(800.0/pl);
88  G4double cpt = std::max(6.23 * std::exp(-1.79*x), 0.3);
89  G4double alphac = 100.0 * 1.0e-6;
90  G4double aaa = (1 + apt) * (1 - std::exp(-btmax))/b;
91  G4double argu = psq * alphac;
92 
93  if(argu >= 8) {
94  argu = 0.0;
95  } else {
96  argu = std::exp(-4.0 * argu);
97  }
98 
99  G4double aac = cpt * (1.0 - argu)/alphac;
100  G4double fracpn = aaa/(aac + aaa);
101  if(Random::shoot() > fracpn) {
102  z = std::exp(-4.0 * psq *alphac);
103  iexpi = 1;
104  y = 1.0 - ranres*(1.0 - z);
105  T = std::log(y)/alphac;
106  }
107  }
108  }
109 
110  G4double ctet = 1.0 + 0.5*T/psq;
111  if(std::abs(ctet) > 1.0) ctet = Math::sign(ctet);
112  G4double stet = std::sqrt(1.0 - ctet*ctet);
113  G4double rndm = Random::shoot();
114 
115  G4double fi = Math::twoPi * rndm;
116  G4double cfi = std::cos(fi);
117  G4double sfi = std::sin(fi);
118 
120  G4double zz = std::pow(particle1->getMomentum().getZ(), 2);
121 
122  if(xx >= (zz * 1.0e-8)) {
124  G4double yn = std::sqrt(xx);
125  G4double zn = yn * pnorm;
126  G4double ex[3], ey[3], ez[3];
127  ez[0] = p.getX() / pnorm;
128  ez[1] = p.getY() / pnorm;
129  ez[2] = p.getZ() / pnorm;
130 
131  // Vector Ex is chosen arbitrarily:
132  ex[0] = p.getY() / yn;
133  ex[1] = -p.getX() / yn;
134  ex[2] = 0.0;
135 
136  ey[0] = p.getX() * p.getZ() / zn;
137  ey[1] = p.getY() * p.getZ() / zn;
138  ey[2] = -xx/zn;
139 
140  G4double pX = (ex[0]*cfi*stet + ey[0]*sfi*stet + ez[0]*ctet) * pnorm;
141  G4double pY = (ex[1]*cfi*stet + ey[1]*sfi*stet + ez[1]*ctet) * pnorm;
142  G4double pZ = (ex[2]*cfi*stet + ey[2]*sfi*stet + ez[2]*ctet) * pnorm;
143 
144  ThreeVector p1momentum = ThreeVector(pX, pY, pZ);
145  particle1->setMomentum(p1momentum);
146  particle2->setMomentum(-p1momentum);
147  } else { // if(xx < (zz * 1.0e-8)) {
148  G4double momZ = particle1->getMomentum().getZ();
149  G4double pX = momZ * cfi * stet;
150  G4double pY = momZ * sfi * stet;
151  G4double pZ = momZ * ctet;
152 
153  ThreeVector p1momentum(pX, pY, pZ);
154  particle1->setMomentum(p1momentum);
155  particle2->setMomentum(-p1momentum);
156  }
157 
158  // Handle backward scattering here.
159 
160  if((particle1->getType() == Proton && particle2->getType() == Neutron) ||
161  (particle1->getType() == Neutron && particle2->getType() == Proton)) {
162  rndm = Random::shoot();
163  apt = 1.0;
164  if(pl > 800.0) {
165  apt = std::pow(800.0/pl, 2);
166  }
167  if(iexpi == 1 || rndm > 1.0/(1.0 + apt)) {
168  particle1->setType(p2TypeOld);
169  particle2->setType(p1TypeOld);
170  }
171  }
172 
173  // Note: there is no need to update the kinetic energies of the particles,
174  // as this is elastic scattering.
175 
178 
179  }
180 
181 }