Geant4-11
Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes
G4INCL::DeltaProductionChannel Class Reference

#include <G4INCLDeltaProductionChannel.hh>

Inheritance diagram for G4INCL::DeltaProductionChannel:
G4INCL::IChannel

Public Member Functions

 DeltaProductionChannel (Particle *, Particle *)
 
void fillFinalState (FinalState *fs)
 
FinalState * getFinalState ()
 
virtual ~DeltaProductionChannel ()
 

Private Member Functions

G4double sampleDeltaMass (G4double ecm)
 

Private Attributes

Particle * particle1
 
Particle * particle2
 

Static Private Attributes

static const G4int maxTries = 100000
 

Detailed Description

Definition at line 48 of file G4INCLDeltaProductionChannel.hh.

Constructor & Destructor Documentation

◆ DeltaProductionChannel()

G4INCL::DeltaProductionChannel::DeltaProductionChannel ( Particle *  p1,
Particle *  p2 
)

Definition at line 49 of file G4INCLDeltaProductionChannel.cc.

51 : particle1(p1), particle2(p2)
52 {}

◆ ~DeltaProductionChannel()

G4INCL::DeltaProductionChannel::~DeltaProductionChannel ( )
virtual

Definition at line 54 of file G4INCLDeltaProductionChannel.cc.

54{}

Member Function Documentation

◆ fillFinalState()

void G4INCL::DeltaProductionChannel::fillFinalState ( FinalState *  fs)
virtual

Delta production

The production is not isotropic in this version it has the same exp(b*t) structure as the nn elastic scattering (formula 2.3 of j.cugnon et al, nucl phys a352(1981)505) parametrization of b taken from ref. prc56(1997)2431

Implements G4INCL::IChannel.

Definition at line 95 of file G4INCLDeltaProductionChannel.cc.

95 {
104 // 100 IF (K4.NE.1) GO TO 101 // ThA K4 = 2 by default
105 // ParticleType p1TypeOld = particle1->getType();
106 // ParticleType p2TypeOld = particle2->getType();
107 G4double ecm = KinematicsUtils::totalEnergyInCM(particle1, particle2);
108
109 const G4int isospin = ParticleTable::getIsospin(particle1->getType()) +
110 ParticleTable::getIsospin(particle2->getType());
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
119 G4double pnorm = KinematicsUtils::momentumInCM(ecm, ParticleTable::effectiveNucleonMass, xmdel);
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;
132 G4double x = 0.001 * KinematicsUtils::momentumInLab(ecm*ecm, ParticleTable::effectiveNucleonMass, 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
180 +ParticleTable::effectiveNucleonMass2);
181 // if(k4.ne.0) go to 161
182
183 // long-lived delta
184 if (index != 1) {
185 ThreeVector mom(xp1, xp2, xp3);
186 particle1->setMomentum(mom);
187 // e1=ecm-eout1
188 } else {
189 ThreeVector mom(-xp1, -xp2, -xp3);
190 particle1->setMomentum(mom);
191 // e1=ecm-eout1
192 }
193
194 particle1->setEnergy(ecm - e3);
195 particle2->setEnergy(e3);
196 particle2->setMomentum(-particle1->getMomentum());
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)
202 G4int is1 = ParticleTable::getIsospin(particle1->getType());
203 G4int is2 = ParticleTable::getIsospin(particle2->getType());
204 if (isospin == 0) {
205 if(index2 == 1) {
206 G4int isi=is1;
207 is1=is2;
208 is2=isi;
209 }
210 particle1->setHelicity(0.0);
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
220 if(is1 == ParticleTable::getIsospin(DeltaMinus)) {
221 particle1->setType(DeltaMinus);
222 } else if(is1 == ParticleTable::getIsospin(DeltaZero)) {
223 particle1->setType(DeltaZero);
224 } else if(is1 == ParticleTable::getIsospin(DeltaPlus)) {
225 particle1->setType(DeltaPlus);
226 } else if(is1 == ParticleTable::getIsospin(DeltaPlusPlus)) {
227 particle1->setType(DeltaPlusPlus);
228 }
229
230 if(is2 == ParticleTable::getIsospin(Proton)) {
231 particle2->setType(Proton);
232 } else if(is2 == ParticleTable::getIsospin(Neutron)) {
233 particle2->setType(Neutron);
234 }
235
236 if(particle1->isDelta()) particle1->setMass(xmdel);
237 if(particle2->isDelta()) particle2->setMass(xmdel);
238
239 fs->addModifiedParticle(particle1);
240 fs->addModifiedParticle(particle2);
241 }
e3
static const G4double e3[45]
Definition: G4BarashenkovData.hh:176
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4INCL::Particle::setMass
void setMass(G4double mass)
Definition: G4INCLParticle.hh:778
G4INCL::Particle::setHelicity
void setHelicity(G4double h)
Definition: G4INCLParticle.hh:827
G4INCL::Particle::getMomentum
const G4INCL::ThreeVector & getMomentum() const
Definition: G4INCLParticle.hh:794
G4INCL::Particle::setMomentum
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
Definition: G4INCLParticle.hh:808
G4INCL::Particle::getType
G4INCL::ParticleType getType() const
Definition: G4INCLParticle.hh:178
G4INCL::Particle::setEnergy
void setEnergy(G4double energy)
Definition: G4INCLParticle.hh:786
G4INCL::Particle::setType
void setType(ParticleType t)
Definition: G4INCLParticle.hh:187
G4INCL::Particle::isDelta
G4bool isDelta() const
Is it a Delta?
Definition: G4INCLParticle.hh:361
G4INCL::ThreeVector::mag
G4double mag() const
Definition: G4INCLThreeVector.hh:73
G4INCL::KinematicsUtils::totalEnergyInCM
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
Definition: G4INCLKinematicsUtils.cc:94
G4INCL::KinematicsUtils::momentumInLab
G4double momentumInLab(Particle const *const p1, Particle const *const p2)
gives the momentum in the lab frame of two particles.
Definition: G4INCLKinematicsUtils.cc:135
G4INCL::KinematicsUtils::momentumInCM
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
Definition: G4INCLKinematicsUtils.cc:107
G4INCL::Math::twoPi
const G4double twoPi
Definition: G4INCLGlobals.hh:69
G4INCL::Math::sign
G4int sign(const T t)
Definition: G4INCLGlobals.hh:107
G4INCL::ParticleTable::effectiveNucleonMass2
const G4double effectiveNucleonMass2
Definition: G4INCLParticleTable.hh:70
G4INCL::ParticleTable::getIsospin
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
Definition: G4INCLParticleTable.cc:478
G4INCL::ParticleTable::effectiveNucleonMass
const G4double effectiveNucleonMass
Definition: G4INCLParticleTable.hh:69
G4INCL::Random::shoot
G4double shoot()
Definition: G4INCLRandom.cc:93

References G4INCL::FinalState::addModifiedParticle(), G4INCL::DeltaMinus, G4INCL::DeltaPlus, G4INCL::DeltaPlusPlus, G4INCL::DeltaZero, e3, G4INCL::ParticleTable::effectiveNucleonMass, G4INCL::ParticleTable::effectiveNucleonMass2, G4INCL::ParticleTable::getIsospin(), G4INCL::Particle::getMomentum(), G4INCL::Particle::getType(), G4INCL::ThreeVector::getX(), G4INCL::ThreeVector::getY(), G4INCL::ThreeVector::getZ(), G4INCL::Particle::isDelta(), G4INCL::ThreeVector::mag(), G4INCL::KinematicsUtils::momentumInCM(), G4INCL::KinematicsUtils::momentumInLab(), G4INCL::Neutron, particle1, particle2, G4INCL::ThreeVector::perp2(), G4INCL::Proton, sampleDeltaMass(), G4INCL::Particle::setEnergy(), G4INCL::Particle::setHelicity(), G4INCL::Particle::setMass(), G4INCL::Particle::setMomentum(), G4INCL::Particle::setType(), G4INCL::Random::shoot(), G4INCL::Math::sign(), G4INCL::KinematicsUtils::totalEnergyInCM(), and G4INCL::Math::twoPi.

◆ getFinalState()

FinalState * G4INCL::IChannel::getFinalState ( )
inherited

Definition at line 50 of file G4INCLIChannel.cc.

50 {
51 FinalState *fs = new FinalState;
52 fillFinalState(fs);
53 return fs;
54 }
G4INCL::IChannel::fillFinalState
virtual void fillFinalState(FinalState *fs)=0

References G4INCL::IChannel::fillFinalState().

◆ sampleDeltaMass()

G4double G4INCL::DeltaProductionChannel::sampleDeltaMass ( G4double  ecm)
private

Definition at line 56 of file G4INCLDeltaProductionChannel.cc.

56 {
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');
74 return ParticleTable::minDeltaMass;
75 }
76
77 G4double rndm = ParticleTable::minDeltaMassRndm + Random::shoot() * deltaMassRndmRange;
78 y = std::tan(rndm);
79 x = ParticleTable::effectiveDeltaMass + 0.5*ParticleTable::effectiveDeltaWidth*y;
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 }
INCL_WARN
#define INCL_WARN(x)
Definition: G4INCLLogger.hh:232
G4bool
bool G4bool
Definition: G4Types.hh:86
G4INCL::ParticleTable::effectiveDeltaWidth
const G4double effectiveDeltaWidth
Definition: G4INCLParticleTable.hh:72
G4INCL::ParticleTable::effectiveDeltaMass
const G4double effectiveDeltaMass
Definition: G4INCLParticleTable.hh:71
G4INCL::ParticleTable::minDeltaMass
G4ThreadLocal G4double minDeltaMass
Definition: G4INCLParticleTable.cc:1232
G4INCL::ParticleTable::minDeltaMassRndm
G4ThreadLocal G4double minDeltaMassRndm
Definition: G4INCLParticleTable.cc:1234

References G4INCL::ParticleTable::effectiveDeltaMass, G4INCL::ParticleTable::effectiveDeltaWidth, G4INCL::ParticleTable::effectiveNucleonMass, INCL_WARN, maxTries, G4INCL::ParticleTable::minDeltaMass, G4INCL::ParticleTable::minDeltaMassRndm, and G4INCL::Random::shoot().

Referenced by fillFinalState().

Field Documentation

◆ maxTries

const G4int G4INCL::DeltaProductionChannel::maxTries = 100000
staticprivate

Definition at line 60 of file G4INCLDeltaProductionChannel.hh.

Referenced by sampleDeltaMass().

◆ particle1

Particle* G4INCL::DeltaProductionChannel::particle1
private

Definition at line 58 of file G4INCLDeltaProductionChannel.hh.

Referenced by fillFinalState().

◆ particle2

Particle * G4INCL::DeltaProductionChannel::particle2
private

Definition at line 58 of file G4INCLDeltaProductionChannel.hh.

Referenced by fillFinalState().

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