Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4INCLPionResonanceDecayChannel.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 
39 #include "G4INCLKinematicsUtils.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
44 // #include <cassert>
45 
46 
47 namespace G4INCL {
48 
49  const G4double PionResonanceDecayChannel::angularSlope = 0.; // Slope to be defined, if needed.
50 
52  :theParticle(p), incidentDirection(dir)
53  { }
54 
56 
57 // Decay during the intranuclear cascade for Omega only
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  }
78 
79  void PionResonanceDecayChannel::sampleAngles(G4double *ctet_par, G4double *stet_par, G4double *phi_par) {
80 
81  (*ctet_par) = -1.0 + 2.0*Random::shoot();
82  if(std::abs(*ctet_par) > 1.0) (*ctet_par) = Math::sign(*ctet_par);
83  (*stet_par) = std::sqrt(1.-(*ctet_par)*(*ctet_par));
84  (*phi_par) = Math::twoPi * Random::shoot();
85  }
86 
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
95  G4double drnd=Random::shoot();
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;
102  createdType = Photon;
103  }
104  else if (drnd < 0.7265) {
105 // 3 pi0
107  pionType1 = PiZero;
108  pionType2 = PiZero;
109  }
110  else if (drnd < 0.9575) {
111 // pi+ pi- pi0
113  pionType1 = PiPlus;
114  pionType2 = PiMinus;
115  }
116  else {
117 // pi+ pi- 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
127  pionType1 = PiPlus;
128  pionType2 = PiMinus;
129  }
130  else if (drnd < 0.9845) {
131 // pi0 photon
132  nbpart=2;
134  createdType = Photon;
135  }
136  else {
137 // pi+ pi-
138  nbpart=2;
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);
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) {
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);
188 
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 
205  fs->addCreatedParticle(Pion1);
206  fs->addCreatedParticle(Pion2);
207 
208  // PhaseSpaceGenerator::generateBiased(sqrtS, list, 0, angularSlope); Biasing?
209  PhaseSpaceGenerator::generate(sqrtS, list);
210  }
211 
212  }
213 }
void addModifiedParticle(Particle *p)
const G4INCL::ThreeVector & getPosition() const
G4double getEnergy() const
TTree * t1
Definition: plottest35.C:26
const char * p
Definition: xmltok.h:285
G4double getMass() const
Get the cached particle mass.
Double_t beta
PionResonanceDecayChannel(Particle *, ThreeVector const &)
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
G4double getZ() const
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
static constexpr double m
Definition: G4SIunits.hh:129
void sampleAngles(G4double *, G4double *, G4double *)
double G4double
Definition: G4Types.hh:76
G4double shoot()
Definition: G4INCLRandom.cc:93
G4double getINCLMass(const G4int A, const G4int Z)
Get INCL nuclear mass (in MeV/c^2)
G4double getX() const
G4double perp() const
void setType(ParticleType t)
G4int sign(const T t)
const G4double twoPi
TTree * t2
Definition: plottest35.C:36
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
int G4int
Definition: G4Types.hh:78
#define INCL_FATAL(x)
TDirectory * dir
G4INCL::ParticleType getType() const
G4double mag() const
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
void addCreatedParticle(Particle *p)
G4double getY() const
const G4double hc
[MeV*fm]