Geant4-11
G4INCLDeltaProductionChannel.cc
Go to the documentation of this file.
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
41#include "G4INCLRandom.hh"
42#include "G4INCLGlobals.hh"
43#include "G4INCLLogger.hh"
44
45namespace G4INCL {
46
48
50 Particle *p2)
51 : particle1(p1), particle2(p2)
52 {}
53
55
57 const G4double maxDeltaMass = ecm - ParticleTable::effectiveNucleonMass - 1.0;
58 const G4double maxDeltaMassRndm = std::atan((maxDeltaMass-ParticleTable::effectiveDeltaMass)*2./ParticleTable::effectiveDeltaWidth);
59 const G4double deltaMassRndmRange = maxDeltaMassRndm - ParticleTable::minDeltaMassRndm;
60// assert(deltaMassRndmRange>0.);
61
62 G4double y=ecm*ecm;
63 G4double q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
64 G4double q3=std::pow(std::sqrt(q2), 3.);
65 const G4double f3max=q3/(q3+5.832E6); // 5.832E6 = 180^3
66 G4double x;
67
68 G4int nTries = 0;
69 G4bool success = false;
70 while(!success) { /* Loop checking, 10.07.2015, D.Mancusi */
71 if(++nTries >= maxTries) {
72 INCL_WARN("DeltaProductionChannel::sampleDeltaMass loop was stopped because maximum number of tries was reached. Minimum delta mass "
73 << ParticleTable::minDeltaMass << " MeV with CM energy " << ecm << " MeV may be unphysical." << '\n');
75 }
76
77 G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange;
78 y = std::tan(rndm);
80// assert(x>=ParticleTable::minDeltaMass && ecm >= x + ParticleTable::effectiveNucleonMass + 1.0);
81
82 // generation of the delta mass with the penetration factor
83 // (see prc56(1997)2431)
84 y=x*x;
85 q2=(y-1.157776E6)*(y-6.4E5)/y/4.0; // 1.157776E6 = 1076^2, 6.4E5 = 800^2
86 q3=std::pow(std::sqrt(q2), 3.);
87 const G4double f3=q3/(q3+5.832E6); // 5.832E6 = 180^3
88 rndm = Random::shoot();
89 if (rndm*f3max < f3)
90 success = true;
91 }
92 return x;
93 }
94
104 // 100 IF (K4.NE.1) GO TO 101 // ThA K4 = 2 by default
105 // ParticleType p1TypeOld = particle1->getType();
106 // ParticleType p2TypeOld = particle2->getType();
108
111
112 // Calculate the outcome of the channel:
113 const ThreeVector &particle1Momentum = particle1->getMomentum();
114 G4double pin = particle1Momentum.mag();
115 G4double rndm = 0.0, b = 0.0;
116
117 G4double xmdel = sampleDeltaMass(ecm);
118 // deltaProduction103: // This label is not used
120 if (pnorm <= 0.0) pnorm=0.000001;
121 G4int index=0;
122 G4int index2=0;
123 rndm = Random::shoot();
124 if (rndm < 0.5) index=1;
125 if (isospin == 0) { // pn case
126 rndm = Random::shoot();
127 if (rndm < 0.5) index2=1;
128 }
129
130 // G4double x=0.001*0.5*ecm*std::sqrt(ecm*ecm-4.*ParticleTable::effectiveNucleonMass2)
131 // / ParticleTable::effectiveNucleonMass;
133 if(x < 1.4) {
134 b=(5.287/(1.+std::exp((1.3-x)/0.05)))*1.e-6;
135 } else {
136 b=(4.65+0.706*(x-1.4))*1.e-6;
137 }
138 G4double xkh = 2.*b*pin*pnorm;
139 rndm = Random::shoot();
140 G4double ctet=1.0+std::log(1.-rndm*(1.-std::exp(-2.*xkh)))/xkh;
141 if(std::abs(ctet) > 1.0) ctet = Math::sign(ctet);
142 G4double stet = std::sqrt(1.-ctet*ctet);
143
144 rndm = Random::shoot();
145 G4double fi = Math::twoPi*rndm;
146 G4double cfi = std::cos(fi);
147 G4double sfi = std::sin(fi);
148 // delta production: correction of the angular distribution 02/09/02
149
150 G4double xx = particle1Momentum.perp2();
151 const G4double particle1MomentumZ = particle1Momentum.getZ();
152 G4double zz = std::pow(particle1MomentumZ, 2);
153 G4double xp1, xp2, xp3;
154 if (xx >= zz*1.e-8) {
155 G4double yn = std::sqrt(xx);
156 G4double zn = yn*pin;
157 G4double ex[3], ey[3], ez[3];
158 G4double p1 = particle1Momentum.getX();
159 G4double p2 = particle1Momentum.getY();
160 G4double p3 = particle1MomentumZ;
161 ez[0] = p1/pin;
162 ez[1] = p2/pin;
163 ez[2] = p3/pin;
164 ex[0] = p2/yn;
165 ex[1] = -p1/yn;
166 ex[2] = 0.0;
167 ey[0] = p1*p3/zn;
168 ey[1] = p2*p3/zn;
169 ey[2] = -xx/zn;
170 xp1 = (ex[0]*cfi*stet+ey[0]*sfi*stet+ez[0]*ctet)*pnorm;
171 xp2 = (ex[1]*cfi*stet+ey[1]*sfi*stet+ez[1]*ctet)*pnorm;
172 xp3 = (ex[2]*cfi*stet+ey[2]*sfi*stet+ez[2]*ctet)*pnorm;
173 }else {
174 xp1=pnorm*stet*cfi;
175 xp2=pnorm*stet*sfi;
176 xp3=pnorm*ctet;
177 }
178 // end of correction angular distribution of delta production
179 G4double e3 = std::sqrt(xp1*xp1+xp2*xp2+xp3*xp3
181 // if(k4.ne.0) go to 161
182
183 // long-lived delta
184 if (index != 1) {
185 ThreeVector mom(xp1, xp2, xp3);
187 // e1=ecm-eout1
188 } else {
189 ThreeVector mom(-xp1, -xp2, -xp3);
191 // e1=ecm-eout1
192 }
193
194 particle1->setEnergy(ecm - e3);
197
198 // SYMMETRIZATION OF CHARGES IN pn -> N DELTA
199 // THE TEST ON "INDEX" ABOVE SYMETRIZES THE EXCITATION OF ONE
200 // OF THE NUCLEONS WITH RESPECT TO THE DELTA EXCITATION
201 // (SEE NOTE 16/10/97)
204 if (isospin == 0) {
205 if(index2 == 1) {
206 G4int isi=is1;
207 is1=is2;
208 is2=isi;
209 }
211 } else {
212 rndm = Random::shoot();
213 if (rndm >= 0.25) {
214 is1=3*is1;
215 is2=-is2;
216 }
217 particle1->setHelicity(ctet*ctet);
218 }
219
222 } else if(is1 == ParticleTable::getIsospin(DeltaZero)) {
224 } else if(is1 == ParticleTable::getIsospin(DeltaPlus)) {
226 } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus)) {
228 }
229
232 } else if(is2 == ParticleTable::getIsospin(Neutron)) {
234 }
235
236 if(particle1->isDelta()) particle1->setMass(xmdel);
237 if(particle2->isDelta()) particle2->setMass(xmdel);
238
241 }
242}
static const G4double e3[45]
#define INCL_WARN(x)
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
void addModifiedParticle(Particle *p)
void setMass(G4double mass)
void setHelicity(G4double h)
const G4INCL::ThreeVector & getMomentum() const
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
G4INCL::ParticleType getType() const
void setEnergy(G4double energy)
void setType(ParticleType t)
G4bool isDelta() const
Is it a Delta?
G4double getY() const
G4double getZ() const
G4double mag() const
G4double perp2() const
G4double getX() const
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
G4double momentumInLab(Particle const *const p1, Particle const *const p2)
gives the momentum in the lab frame of two particles.
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
const G4double twoPi
G4int sign(const T t)
const G4double effectiveDeltaWidth
const G4double effectiveDeltaMass
G4ThreadLocal G4double minDeltaMass
G4ThreadLocal G4double minDeltaMassRndm
const G4double effectiveNucleonMass2
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
const G4double effectiveNucleonMass
G4double shoot()
Definition: G4INCLRandom.cc:93