Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCL::PionResonanceDecayChannel Class Reference

#include <G4INCLPionResonanceDecayChannel.hh>

+ Inheritance diagram for G4INCL::PionResonanceDecayChannel:

Public Member Functions

 PionResonanceDecayChannel (Particle *, ThreeVector const &)
 
virtual ~PionResonanceDecayChannel ()
 
void fillFinalState (FinalState *fs)
 
- Public Member Functions inherited from G4INCL::IChannel
 IChannel ()
 
virtual ~IChannel ()
 
FinalStategetFinalState ()
 

Static Public Member Functions

static G4double computeDecayTime (Particle *p)
 

Detailed Description

Definition at line 49 of file G4INCLPionResonanceDecayChannel.hh.

Constructor & Destructor Documentation

◆ PionResonanceDecayChannel()

G4INCL::PionResonanceDecayChannel::PionResonanceDecayChannel ( Particle * p,
ThreeVector const & dir )

Definition at line 51 of file G4INCLPionResonanceDecayChannel.cc.

52 :theParticle(p), incidentDirection(dir)
53 { }

◆ ~PionResonanceDecayChannel()

G4INCL::PionResonanceDecayChannel::~PionResonanceDecayChannel ( )
virtual

Definition at line 55 of file G4INCLPionResonanceDecayChannel.cc.

55{}

Member Function Documentation

◆ computeDecayTime()

G4double G4INCL::PionResonanceDecayChannel::computeDecayTime ( Particle * p)
static

Definition at line 58 of file G4INCLPionResonanceDecayChannel.cc.

58 {
59 const G4double m = p->getMass();
60 const G4double geff = p->getEnergy()/m;
61// const G4double geta = 1.31e-3;
62 const G4double gomega = 8.49;
63 G4double gg=0.;
64 switch (p->getType()) {
65/* case Eta:
66 gg=geta;
67 break;*/
68 case Omega:
69 gg=gomega;
70 break;
71 default:
72 INCL_FATAL("Unrecognized pion resonance type; type=" << p->getType() << '\n');
73 break;
74 }
75 const G4double tpires = -G4INCL::PhysicalConstants::hc/gg*std::log(Random::shoot())*geff;
76 return tpires;
77 }
#define INCL_FATAL(x)
double G4double
Definition G4Types.hh:83
const G4double hc
[MeV*fm]
G4double shoot()

Referenced by G4INCL::StandardPropagationModel::generateDecays().

◆ fillFinalState()

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

Implements G4INCL::IChannel.

Definition at line 87 of file G4INCLPionResonanceDecayChannel.cc.

87 {
88
89 ParticleType createdType;
90 ParticleType pionType1=Neutron; // to avoid forgetting pionType definition when 3 particles are emitted
91 ParticleType pionType2=Neutron;
92
93 const G4double sqrtS = theParticle->getMass();
94 G4int nbpart = 3; // number of emitted particles
96 switch (theParticle->getType()) {
97 case Eta:
98 if (drnd < 0.3972) { // renormalized to the only four decays taken into account here
99// 2 photons
100 nbpart=2;
101 theParticle->setType(Photon);
102 createdType = Photon;
103 }
104 else if (drnd < 0.7265) {
105// 3 pi0
106 theParticle->setType(PiZero);
107 pionType1 = PiZero;
108 pionType2 = PiZero;
109 }
110 else if (drnd < 0.9575) {
111// pi+ pi- pi0
112 theParticle->setType(PiZero);
113 pionType1 = PiPlus;
114 pionType2 = PiMinus;
115 }
116 else {
117// pi+ pi- photon
118 theParticle->setType(Photon);
119 pionType1 = PiPlus;
120 pionType2 = PiMinus;
121 }
122 break;
123 case Omega:
124 if (drnd < 0.9009) { // renormalized to the only three decays taken into account here
125// pi+ pi- pi0
126 theParticle->setType(PiZero);
127 pionType1 = PiPlus;
128 pionType2 = PiMinus;
129 }
130 else if (drnd < 0.9845) {
131// pi0 photon
132 nbpart=2;
133 theParticle->setType(PiZero);
134 createdType = Photon;
135 }
136 else {
137// pi+ pi-
138 nbpart=2;
139 theParticle->setType(PiPlus);
140 createdType = PiMinus;
141 }
142 break;
143 default:
144 INCL_FATAL("Unrecognized pion resonance type; type=" << theParticle->getType() << '\n');
145 break;
146 }
147
148 if (nbpart == 2) {
149 G4double fi, ctet, stet;
150 sampleAngles(&ctet, &stet, &fi);
151
152 G4double cfi = std::cos(fi);
153 G4double sfi = std::sin(fi);
154 G4double beta = incidentDirection.mag();
155
156 G4double q1, q2, q3;
157 G4double sal=0.0;
158 if (beta >= 1.0e-10)
159 sal = incidentDirection.perp()/beta;
160 if (sal >= 1.0e-6) {
161 G4double b1 = incidentDirection.getX();
162 G4double b2 = incidentDirection.getY();
163 G4double b3 = incidentDirection.getZ();
164 G4double cal = b3/beta;
165 G4double t1 = ctet+cal*stet*sfi/sal;
166 G4double t2 = stet/sal;
167 q1=(b1*t1+b2*t2*cfi)/beta;
168 q2=(b2*t1-b1*t2*cfi)/beta;
169 q3=(b3*t1/beta-t2*sfi);
170 } else {
171 q1 = stet*cfi;
172 q2 = stet*sfi;
173 q3 = ctet;
174 }
175
177 theParticle->getMass(),
178 ParticleTable::getINCLMass(createdType));
179 q1 *= xq;
180 q2 *= xq;
181 q3 *= xq;
182
183 ThreeVector createdMomentum(q1, q2, q3);
184 ThreeVector createdPosition(theParticle->getPosition());
185 Particle *createdParticle = new Particle(createdType, createdMomentum, createdPosition);
186 theParticle->setMomentum(-createdMomentum);
187 theParticle->adjustEnergyFromMomentum();
188
189 fs->addModifiedParticle(theParticle);
190 fs->addCreatedParticle(createdParticle);
191
192 }
193 else if (nbpart == 3) {
194// assert(pionType1!=Neutron && pionType2!=Neutron);
195 ParticleList list;
196 list.push_back(theParticle);
197 const ThreeVector &rposdecay = theParticle->getPosition();
198 const ThreeVector zero;
199 Particle *Pion1 = new Particle(pionType1,zero,rposdecay);
200 Particle *Pion2 = new Particle(pionType2,zero,rposdecay);
201 list.push_back(Pion1);
202 list.push_back(Pion2);
203
204 fs->addModifiedParticle(theParticle);
205 fs->addCreatedParticle(Pion1);
206 fs->addCreatedParticle(Pion2);
207
208 // PhaseSpaceGenerator::generateBiased(sqrtS, list, 0, angularSlope); Biasing?
210 }
211
212 }
int G4int
Definition G4Types.hh:85
const G4INCL::ThreeVector & getPosition() const
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
G4INCL::ParticleType getType() const
void setType(ParticleType t)
G4double getMass() const
Get the cached particle mass.
G4double getY() const
G4double getZ() const
G4double perp() const
G4double getX() const
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
G4double getINCLMass(const G4int A, const G4int Z, const G4int S)
Get INCL nuclear mass (in MeV/c^2)
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.

The documentation for this class was generated from the following files: