Geant4-11
Public Member Functions | Static Private Attributes
G4NRESP71M03 Class Reference

#include <G4NRESP71M03.hh>

Public Member Functions

G4int ApplyMechanismABE (G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds)
 
G4int ApplyMechanismI_NBeA2A (G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI)
 
G4int ApplyMechanismII_ACN2A (G4ReactionProduct &neut, G4ReactionProduct &carb, G4ReactionProduct *theProds, const G4double QI)
 
void DKINMA (G4ReactionProduct *p1, G4ReactionProduct *p2, G4ReactionProduct *p3, G4ReactionProduct *p4, const G4double Q, const G4double costhcm3)
 
 G4NRESP71M03 ()
 
 ~G4NRESP71M03 ()
 

Static Private Attributes

static const G4double B2 [ndist][nrhos]
 
static const G4double BEN2 [ndist] = { 5700., 8000., 8640., 8990., 9220., 9410., 9830., 10400., 10800., 11250., 11460., 11870., 12140., 12320., 12570., 12940., 13420., 13760., 14020., 14200., 14440., 14620., 14820., 15050., 15660., 15980., 16470., 16940., 17970., 18000., 19000., 20000. }
 
static const G4int ndist = 32
 
static const G4int nrhos = 51
 

Detailed Description

Definition at line 36 of file G4NRESP71M03.hh.

Constructor & Destructor Documentation

◆ G4NRESP71M03()

G4NRESP71M03::G4NRESP71M03 ( )
inline

Definition at line 40 of file G4NRESP71M03.hh.

40{;};

◆ ~G4NRESP71M03()

G4NRESP71M03::~G4NRESP71M03 ( )
inline

Definition at line 41 of file G4NRESP71M03.hh.

41{;};

Member Function Documentation

◆ ApplyMechanismABE()

G4int G4NRESP71M03::ApplyMechanismABE ( G4ReactionProduct &  neut,
G4ReactionProduct &  carb,
G4ReactionProduct *  theProds 
)

Definition at line 202 of file G4NRESP71M03.cc.

203 {
204 G4double CosThetaCMAlpha = 0.; // Cosine of the angle of emission of the alpha particle (theta).
205
206 G4double Kn = neut.GetKineticEnergy(); // Neutron energy in the center of mass system.
207 if ( Kn > 5.7*MeV )
208 {
209 // Sorting.
210 for ( G4int i=1; i<ndist; i++ )
211 {
212 if ( BEN2[i] >= Kn/keV )
213 {
214 // Ok. The neutron energy is between values i-1 and i of BEN2. Taking them.
215 G4double BE1 = BEN2[i-1];
216 G4double BE2 = BEN2[i];
217
218 // Performing energy and angle interpolation.
219
220 G4double FRA = G4UniformRand()*49.99999999; // Sorting the bin of the cumulative probability FRA (Rho). There are 51 values of Rho from 0 to 1 (0; 0.02; 0.04 ... 1).
221 G4double DJTETA = FRA-G4int(FRA); // Distance in bin units (DJTETA) from the low edge of the bin with Rho.
222 G4int JTETA = G4int(FRA)+1; // Getting the bin (JTETA) next to the bin with the value of Rho.
223
224 // Calculating the angle from the cumulative distribution at energy i.
225
226 G4double B11 = B2[i-1][JTETA-1];
227 G4double B12 = B2[i-1][JTETA];
228
229 G4double TCM1 = B11+(B12-B11)*DJTETA; // Angle at energy i.
230
231 // Calculating the angle from the cumulative distribution at energy i-1.
232
233 G4double B21 = B2[i][JTETA-1];
234 G4double B22 = B2[i][JTETA];
235
236 G4double TCM2 = B21+(B22-B21)*DJTETA; // Angle at energy i-1.
237
238 // Interpolating the angle between energy values i and i-1.
239 G4double TCM = (TCM1+(TCM2-TCM1)*(Kn/keV-BE1)/(BE2-BE1))*1.E-4;
240 CosThetaCMAlpha = std::cos(TCM);
241
242 break;
243 }
244 }
245 }
246 else
247 {
248 // Isotropic distribution in CM.
249 CosThetaCMAlpha = 1.-2.*G4UniformRand();
250 }
251
252 //N+12C --> A+9BE
253 theProds[0].SetDefinition(G4Alpha::Alpha());
254 theProds[1].SetDefinition(G4IonTable::GetIonTable()->GetIon(4, 9, 0. ));
255
256 DKINMA(&neut, &carb, &(theProds[0]), &(theProds[1]), -5.71*MeV, CosThetaCMAlpha);
257
258 return 0;
259 }
keV
static constexpr double keV
Definition: G4SIunits.hh:202
MeV
static constexpr double MeV
Definition: G4SIunits.hh:200
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition: G4Alpha.cc:88
G4IonTable::GetIonTable
static G4IonTable * GetIonTable()
Definition: G4IonTable.cc:170
G4NRESP71M03::DKINMA
void DKINMA(G4ReactionProduct *p1, G4ReactionProduct *p2, G4ReactionProduct *p3, G4ReactionProduct *p4, const G4double Q, const G4double costhcm3)
Definition: G4NRESP71M03.cc:76
G4NRESP71M03::B2
static const G4double B2[ndist][nrhos]
Definition: G4NRESP71M03.hh:57
G4NRESP71M03::ndist
static const G4int ndist
Definition: G4NRESP71M03.hh:52
G4NRESP71M03::BEN2
static const G4double BEN2[ndist]
Definition: G4NRESP71M03.hh:55
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:163

References G4Alpha::Alpha(), B2, BEN2, DKINMA(), G4UniformRand, G4IonTable::GetIonTable(), G4ReactionProduct::GetKineticEnergy(), keV, MeV, ndist, and G4ReactionProduct::SetDefinition().

Referenced by G4ParticleHPInelasticCompFS::use_nresp71_model().

◆ ApplyMechanismI_NBeA2A()

G4int G4NRESP71M03::ApplyMechanismI_NBeA2A ( G4ReactionProduct &  neut,
G4ReactionProduct &  carb,
G4ReactionProduct *  theProds,
const G4double  QI 
)

Definition at line 147 of file G4NRESP71M03.cc.

148 {
149 // N+12C --> A+9BE*
150 G4ReactionProduct p4;
151
152 theProds[0].SetDefinition(G4Alpha::Alpha());
153
154 DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2.*G4UniformRand()-1.);
155
156 // 9BE* --> N+8BE
157 G4ReactionProduct p1(p4);
158
159 theProds[1].SetDefinition(G4Neutron::Neutron());
160
161 DKINMA(&p1, NULL, &(theProds[1]), &p4, -QI-7.369, 2.*G4UniformRand()-1.);
162
163 // 8BE --> 2*A
164 p1 = p4;
165
166 theProds[2].SetDefinition(G4Alpha::Alpha());
167 theProds[3].SetDefinition(G4Alpha::Alpha());
168
169 DKINMA(&p1, NULL, &(theProds[2]), &(theProds[3]), 0.09538798439007223351, 2.*G4UniformRand()-1.);
170
171 return 0;
172 }
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103

References G4Alpha::Alpha(), DKINMA(), G4UniformRand, G4Neutron::Neutron(), and G4ReactionProduct::SetDefinition().

Referenced by G4ParticleHPInelasticCompFS::use_nresp71_model().

◆ ApplyMechanismII_ACN2A()

G4int G4NRESP71M03::ApplyMechanismII_ACN2A ( G4ReactionProduct &  neut,
G4ReactionProduct &  carb,
G4ReactionProduct *  theProds,
const G4double  QI 
)

Definition at line 175 of file G4NRESP71M03.cc.

176 {
177 // 12C(N,N')12C'
178 G4ReactionProduct p4;
179
180 theProds[0].SetDefinition(G4Neutron::Neutron());
181
182 DKINMA(&neut, &carb, &(theProds[0]), &p4, QI, 2.*G4UniformRand()-1.);
183
184 // 12C' --> ALPHA+8BE'
185 G4ReactionProduct p1(p4);
186
187 theProds[1].SetDefinition(G4Alpha::Alpha());
188
189 DKINMA(&p1, NULL, &(theProds[1]), &p4, -QI-7.369, 2.*G4UniformRand()-1.);
190
191 // 8BE --> 2*ALPHA
192 p1 = p4;
193
194 theProds[2].SetDefinition(G4Alpha::Alpha());
195 theProds[3].SetDefinition(G4Alpha::Alpha());
196
197 DKINMA(&p1, NULL, &(theProds[2]), &(theProds[3]), 0.09538798439007223351, 2.*G4UniformRand()-1.);
198
199 return 0;
200 }

References G4Alpha::Alpha(), DKINMA(), G4UniformRand, G4Neutron::Neutron(), and G4ReactionProduct::SetDefinition().

Referenced by G4ParticleHPInelasticCompFS::use_nresp71_model().

◆ DKINMA()

void G4NRESP71M03::DKINMA ( G4ReactionProduct *  p1,
G4ReactionProduct *  p2,
G4ReactionProduct *  p3,
G4ReactionProduct *  p4,
const G4double  Q,
const G4double  costhcm3 
)

Definition at line 76 of file G4NRESP71M03.cc.

77 {
78 G4ReactionProduct CM;
79 G4double TotalEnergyCM;
80
81 if ( p2 ) // If it is not a decay reaction...
82 {
83 // Calculating (total momentum, energy and mass) of the center of mass.
84 const G4ThreeVector TotalMomentumLAB = p1->GetMomentum()+p2->GetMomentum();
85 CM.SetMomentum(TotalMomentumLAB);
86
87 const G4double TotalEnergyLAB = p1->GetTotalEnergy()+p2->GetTotalEnergy();
88 CM.SetTotalEnergy(TotalEnergyLAB);
89
90 CM.SetMass(std::sqrt(TotalEnergyLAB*TotalEnergyLAB-TotalMomentumLAB*TotalMomentumLAB));
91
92 // Transforming primary particles from laboratory to center of mass system.
93 p1->Lorentz(*p1, CM);
94 p2->Lorentz(*p2, CM);
95
96 TotalEnergyCM = p1->GetTotalEnergy()+p2->GetTotalEnergy();
97
98 const G4double m4 = (p1->GetMass()+p2->GetMass())-(p3->GetMass()+Q); // Mass of the residual nucleus in the excited state (not in the ground state).
99 p4->SetMass(m4);
100 }
101 else // If it is a decay reaction...
102 {
103 const G4ThreeVector TotalMomentumLAB = p1->GetMomentum();
104 CM.SetMomentum(TotalMomentumLAB);
105
106 const G4double TotalEnergyLAB = p1->GetTotalEnergy();
107 CM.SetTotalEnergy(TotalEnergyLAB);
108
109 CM.SetMass(std::sqrt(TotalEnergyLAB*TotalEnergyLAB-TotalMomentumLAB*TotalMomentumLAB));
110
111 // Transforming primary particles from laboratory to center of mass system (not really necessary in this case).
112 p1->Lorentz(*p1, CM);
113
114 const G4double m4 = p1->GetMass()-(p3->GetMass()+Q); // Mass of the residual nucleus in the excited state (not in the ground state).
115 p4->SetMass(m4);
116
117 TotalEnergyCM = p1->GetTotalEnergy();
118 }
119
120 // Calculating momentum and total energy of the reaction products.
121
122 const G4ThreeVector p1unit = p1->GetMomentum().unit();
123
124 G4RotationMatrix rot(std::acos(p1unit*G4ThreeVector(0., 1., 0.)), std::acos(p1unit*G4ThreeVector(0., 0., 1.)), 0.);
125 rot.inverse();
126
127 const G4double theta = std::acos(costhcm3);
128 const G4double phi = twopi*G4UniformRand();
129
130 const G4double Energy3CM = (std::pow(TotalEnergyCM, 2.)+std::pow(p3->GetMass(), 2.)-std::pow(p4->GetMass(), 2.))/(2.*TotalEnergyCM);
131 p3->SetTotalEnergy(Energy3CM);
132
133 const G4double Momentum3CM = std::sqrt(std::pow(Energy3CM, 2.)-std::pow(p3->GetMass(), 2.));
134 p3->SetMomentum(rot*G4ThreeVector(Momentum3CM*std::sin(theta)*std::cos(phi), Momentum3CM*std::sin(theta)*std::sin(phi), Momentum3CM*costhcm3));
135
136 const G4double Energy4CM = TotalEnergyCM-Energy3CM;
137 p4->SetTotalEnergy(Energy4CM);
138
139 const G4double Momentum4CM = std::sqrt(std::pow(Energy4CM, 2.)-std::pow(p4->GetMass(), 2.));
140 p4->SetMomentum(-Momentum4CM*p3->GetMomentum().unit());
141
142 // Transforming reaction products from center of mass to laboratory system.
143 p3->Lorentz(*p3, -1.*CM);
144 p4->Lorentz(*p4, -1.*CM);
145 }
twopi
static constexpr double twopi
Definition: G4SIunits.hh:56
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:36
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:174
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:141
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:129
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:123
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:209
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:150
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:147
Q
static double Q[]
Definition: nf_gammaFunctions.cc:91

References G4UniformRand, G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ReactionProduct::GetTotalEnergy(), CLHEP::HepRotation::inverse(), G4ReactionProduct::Lorentz(), Q, G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTotalEnergy(), twopi, and CLHEP::Hep3Vector::unit().

Referenced by ApplyMechanismABE(), ApplyMechanismI_NBeA2A(), and ApplyMechanismII_ACN2A().

Field Documentation

◆ B2

const G4double G4NRESP71M03::B2
staticprivate

Definition at line 57 of file G4NRESP71M03.hh.

Referenced by ApplyMechanismABE().

◆ BEN2

const G4double G4NRESP71M03::BEN2 = { 5700., 8000., 8640., 8990., 9220., 9410., 9830., 10400., 10800., 11250., 11460., 11870., 12140., 12320., 12570., 12940., 13420., 13760., 14020., 14200., 14440., 14620., 14820., 15050., 15660., 15980., 16470., 16940., 17970., 18000., 19000., 20000. }
staticprivate

Definition at line 55 of file G4NRESP71M03.hh.

Referenced by ApplyMechanismABE().

◆ ndist

const G4int G4NRESP71M03::ndist = 32
staticprivate

Definition at line 52 of file G4NRESP71M03.hh.

Referenced by ApplyMechanismABE().

◆ nrhos

const G4int G4NRESP71M03::nrhos = 51
staticprivate

Definition at line 53 of file G4NRESP71M03.hh.

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.3