| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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
38 #include "G4INCLPionResonanceDecayChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLPhaseSpaceGenerator.hh"
44 // #include <cassert>
45
46
47 namespace G4INCL {
48
49 const G4double PionResonanceDecayChannel::an
50
51 PionResonanceDecayChannel::PionResonanceDeca
52 :theParticle(p), incidentDirection(dir)
53 { }
54
55 PionResonanceDecayChannel::~PionResonanceDec
56
57 // Decay during the intranuclear cascade for
58 G4double PionResonanceDecayChannel::computeD
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 resonanc
73 break;
74 }
75 const G4double tpires = -G4INCL::PhysicalC
76 return tpires;
77 }
78
79 void PionResonanceDecayChannel::sampleAngles
80
81 (*ctet_par) = -1.0 + 2.0*Random::shoot();
82 if(std::abs(*ctet_par) > 1.0) (*ctet_par)
83 (*stet_par) = std::sqrt(1.-(*ctet_par)*(*c
84 (*phi_par) = Math::twoPi * Random::shoot()
85 }
86
87 void PionResonanceDecayChannel::fillFinalSta
88
89 ParticleType createdType;
90 ParticleType pionType1=Neutron; // to avoi
91 ParticleType pionType2=Neutron;
92
93 const G4double sqrtS = theParticle->getMas
94 G4int nbpart = 3; // number of emitted par
95 G4double drnd=Random::shoot();
96 switch (theParticle->getType()) {
97 case Eta:
98 if (drnd < 0.3972) { // renormalized t
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 t
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 resonanc
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
176 G4double xq = KinematicsUtils::momentumI
177 theParticle->getMa
178 ParticleTable::get
179 q1 *= xq;
180 q2 *= xq;
181 q3 *= xq;
182
183 ThreeVector createdMomentum(q1, q2, q3);
184 ThreeVector createdPosition(theParticle-
185 Particle *createdParticle = new Particle
186 theParticle->setMomentum(-createdMomentu
187 theParticle->adjustEnergyFromMomentum();
188
189 fs->addModifiedParticle(theParticle);
190 fs->addCreatedParticle(createdParticle);
191
192 }
193 else if (nbpart == 3) {
194 // assert(pionType1!=Neutron && pionType2!=Neu
195 ParticleList list;
196 list.push_back(theParticle);
197 const ThreeVector &rposdecay = thePartic
198 const ThreeVector zero;
199 Particle *Pion1 = new Particle(pionType1
200 Particle *Pion2 = new Particle(pionType2
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::generateBia
209 PhaseSpaceGenerator::generate(sqrtS, lis
210 }
211
212 }
213 }
214