| 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 "G4INCLNNbarToNNbar3piChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 #include <algorithm>
45 #include "G4INCLPhaseSpaceGenerator.hh"
46
47 namespace G4INCL {
48
49 NNbarToNNbar3piChannel::NNbarToNNbar3piChann
50 : particle1(p1), particle2(p2)
51 {}
52
53 NNbarToNNbar3piChannel::~NNbarToNNbar3piChan
54
55 void NNbarToNNbar3piChannel::fillFinalState(
56
57 //brief ppbar
58 // p pbar -> p pbar pi+ pi- pi0 (BFMM
59 // p pbar -> p nbar 2pi- pi+ (BFMM 169
60 // p pbar -> n pbar 2pi+ pi- (BFMM 201
61 // p pbar -> n nbar pi+ pi- pi0 (BFMM
62 //
63 //brief npbar
64 // n pbar -> p pbar 2pi- pi+ (same as
65 // n pbar -> p nbar 2pi- pi0 (same as
66 // n pbar -> n nbar 2pi- pi+ (same as
67 // n pbar -> n pbar pi+ pi- pi0 (same
68 //
69 //brief nnbar
70 // n nbar -> n nbar pi+ pi- pi0 (same
71 // n nbar -> p nbar 2pi- pi+ (same as
72 // n nbar -> n pbar 2pi+ pi- (same as
73 // n nbar -> p pbar pi+ pi- pi0 (same
74 //
75 //brief pnbar
76 // p nbar -> p pbar 2pi+ pi- (same as
77 // p nbar -> n pbar 2pi+ pi0 (same as
78 // p nbar -> n nbar 2pi+ pi- (same as
79 // p nbar -> p nbar pi+ pi- pi0 (same
80
81 Particle *nucleon;
82 Particle *antinucleon;
83
84 if(particle1->isNucleon()){
85 nucleon = particle1;
86 antinucleon = particle2;
87 }
88 else{
89 nucleon = particle2;
90 antinucleon = particle1;
91 }
92
93 const G4double plab = 0.001*KinematicsUtil
94 const G4double sqrtS = KinematicsUtils::to
95 const G4double rdm = Random::shoot();
96
97 const std::vector<G4double> BFMM161 =
98 //const G4double Eth_PPbar_PPbar_pip_p
99 const std::vector<G4double> BFMM169 =
100 //const G4double Eth_PPbar_PNbar_2pim_
101 const std::vector<G4double> BFMM201 =
102 //const G4double Eth_PPbar_NPbar_2pip_
103 const std::vector<G4double> BFMM197 =
104 //const G4double Eth_PPbar_NNbar_pip_p
105
106 // pnbar total is same as for npbar
107 // ppbar total is same as for nnbar
108 const G4double totalppbar = KinematicsUtil
109 +KinematicsUtils::compute_xs(BFMM169, plab
110 +KinematicsUtils::compute_xs(BFMM201, plab
111 +KinematicsUtils::compute_xs(BFMM197, plab
112 const G4double totalpnbar = KinematicsUtil
113 +KinematicsUtils::compute_xs(BFMM197, plab
114 +2*KinematicsUtils::compute_xs(BFMM169, pl
115
116 //totalnnbar == totalppbar;
117 //totalpnbar == totalnpbar;
118 ParticleType Pion1;
119 ParticleType Pion2;
120 ParticleType Pion3;
121
122 //setting types of new particles
123 if(nucleon->getType()==Proton){
124 if(antinucleon->getType()==antiProton){
125 if(rdm*totalppbar < KinematicsUtils::c
126 Pion1 = PiMinus;
127 Pion2 = PiPlus;
128 Pion3 = PiZero;
129 if(rdm<0.5){
130 nucleon->setType(Proton);
131 antinucleon->setType(antiProton);
132 }
133 else{
134 nucleon->setType(antiProton);
135 antinucleon->setType(Proton);
136 }
137 }
138 else if(rdm*totalppbar < KinematicsUti
139 +KinematicsUtils::compute_xs(BFMM169,
140 Pion1 = PiMinus;
141 Pion2 = PiMinus;
142 Pion3 = PiPlus;
143 if(rdm<0.5){
144 nucleon->setType(Proton);
145 antinucleon->setType(antiNeutron);
146 }
147 else{
148 nucleon->setType(antiNeutron);
149 antinucleon->setType(Proton);
150 }
151 }
152 else if(rdm*totalppbar < KinematicsUti
153 +KinematicsUtils::compute_xs(std::move
154 +KinematicsUtils::compute_xs(std::move
155 Pion1 = PiPlus;
156 Pion2 = PiPlus;
157 Pion3 = PiMinus;
158 if(rdm<0.5){
159 nucleon->setType(Neutron);
160 antinucleon->setType(antiProton);
161 }
162 else{
163 nucleon->setType(antiProton);
164 antinucleon->setType(Neutron);
165 }
166 }
167 else{ // n nbar pi+ pi- pi0 case
168 Pion1 = PiMinus;
169 Pion2 = PiPlus;
170 Pion3 = PiZero;
171 if(rdm<0.5){
172 nucleon->setType(Neutron);
173 antinucleon->setType(antiNeutron);
174 }
175 else{
176 nucleon->setType(antiNeutron);
177 antinucleon->setType(Neutron);
178 }
179 }
180 }
181 else{ //antiNeutron (pnbar case)
182 if(rdm*totalpnbar < KinematicsUtils::c
183 Pion1 = PiPlus;
184 Pion2 = PiPlus;
185 Pion3 = PiMinus;
186 if(rdm<0.5){
187 nucleon->setType(Proton);
188 antinucleon->setType(antiProton);
189 }
190 else{
191 nucleon->setType(antiProton);
192 antinucleon->setType(Proton);
193 }
194 }
195 else if(rdm*totalppbar < KinematicsUti
196 +KinematicsUtils::compute_xs(BFMM197,
197 Pion1 = PiPlus;
198 Pion2 = PiPlus;
199 Pion3 = PiZero;
200 if(rdm<0.5){
201 nucleon->setType(Neutron);
202 antinucleon->setType(antiProton);
203 }
204 else{
205 nucleon->setType(antiProton);
206 antinucleon->setType(Neutron);
207 }
208 }
209 else if(rdm*totalppbar < 2*KinematicsU
210 +KinematicsUtils::compute_xs(std::move
211 Pion1 = PiPlus;
212 Pion2 = PiPlus;
213 Pion3 = PiMinus;
214 if(rdm<0.5){
215 nucleon->setType(Neutron);
216 antinucleon->setType(antiNeutron);
217 }
218 else{
219 nucleon->setType(antiNeutron);
220 antinucleon->setType(Neutron);
221 }
222 }
223 else{ // p nbar pi+ pi- pi0 case
224 Pion1 = PiMinus;
225 Pion2 = PiPlus;
226 Pion3 = PiZero;
227 if(rdm<0.5){
228 nucleon->setType(Proton);
229 antinucleon->setType(antiNeutron);
230 }
231 else{
232 nucleon->setType(antiNeutron);
233 antinucleon->setType(Proton);
234 }
235 }
236 }
237 }
238 else{ // neutron
239 if(antinucleon->getType()==antiProton){
240 if(rdm*totalpnbar < KinematicsUtils::c
241 Pion1 = PiPlus;
242 Pion2 = PiMinus;
243 Pion3 = PiMinus;
244 if(rdm<0.5){
245 nucleon->setType(Proton);
246 antinucleon->setType(antiProton);
247 }
248 else{
249 nucleon->setType(antiProton);
250 antinucleon->setType(Proton);
251 }
252 }
253 else if(rdm*totalppbar < KinematicsUti
254 +KinematicsUtils::compute_xs(BFMM197,
255 Pion1 = PiMinus;
256 Pion2 = PiMinus;
257 Pion3 = PiZero;
258 if(rdm<0.5){
259 nucleon->setType(Proton);
260 antinucleon->setType(antiNeutron);
261 }
262 else{
263 nucleon->setType(antiNeutron);
264 antinucleon->setType(Proton);
265 }
266 }
267 else if(rdm*totalppbar < 2*KinematicsU
268 +KinematicsUtils::compute_xs(std::move
269 Pion1 = PiPlus;
270 Pion2 = PiMinus;
271 Pion3 = PiMinus;
272 if(rdm<0.5){
273 nucleon->setType(Neutron);
274 antinucleon->setType(antiNeutron);
275 }
276 else{
277 nucleon->setType(antiNeutron);
278 antinucleon->setType(Neutron);
279 }
280 }
281 else{ // n pbar pi+ pi- pi0 case
282 Pion1 = PiMinus;
283 Pion2 = PiPlus;
284 Pion3 = PiZero;
285 if(rdm<0.5){
286 nucleon->setType(Neutron);
287 antinucleon->setType(antiProton);
288 }
289 else{
290 nucleon->setType(antiProton);
291 antinucleon->setType(Neutron);
292 }
293 }
294 }
295 else{ //antiNeutron (nnbar case)
296 if(rdm*totalppbar < KinematicsUtils::c
297 Pion1 = PiMinus;
298 Pion2 = PiPlus;
299 Pion3 = PiZero;
300 if(rdm<0.5){
301 nucleon->setType(Neutron);
302 antinucleon->setType(antiNeutron);
303 }
304 else{
305 nucleon->setType(antiNeutron);
306 antinucleon->setType(Neutron);
307 }
308 }
309 else if(rdm*totalppbar < KinematicsUti
310 +KinematicsUtils::compute_xs(BFMM169,
311 Pion1 = PiMinus;
312 Pion2 = PiMinus;
313 Pion3 = PiPlus;
314 if(rdm<0.5){
315 nucleon->setType(Proton);
316 antinucleon->setType(antiNeutron);
317 }
318 else{
319 nucleon->setType(antiNeutron);
320 antinucleon->setType(Proton);
321 }
322 }
323 else if(rdm*totalppbar < KinematicsUti
324 +KinematicsUtils::compute_xs(std::move
325 +KinematicsUtils::compute_xs(std::move
326 Pion1 = PiPlus;
327 Pion2 = PiPlus;
328 Pion3 = PiMinus;
329 if(rdm<0.5){
330 nucleon->setType(Neutron);
331 antinucleon->setType(antiProton);
332 }
333 else{
334 nucleon->setType(antiProton);
335 antinucleon->setType(Neutron);
336 }
337 }
338 else{ // p pbar pi+ pi- pi0 case
339 Pion1 = PiMinus;
340 Pion2 = PiPlus;
341 Pion3 = PiZero;
342 if(rdm<0.5){
343 nucleon->setType(Proton);
344 antinucleon->setType(antiProton);
345 }
346 else{
347 nucleon->setType(antiProton);
348 antinucleon->setType(Proton);
349 }
350 }
351 }
352 }
353
354 ParticleList list;
355 list.push_back(nucleon);
356 list.push_back(antinucleon);
357 const ThreeVector &rcol = nucleon->getPosi
358 const ThreeVector zero;
359
360 // Create three particle pointers
361 Particle *pion1 = nullptr;
362 Particle *pion2 = nullptr;
363 Particle *pion3 = nullptr;
364
365 // Determine the types of particles based
366 if (rdm < 1.0 / 3.0) {
367 pion1 = new Particle(Pion1, zero, rcol
368 pion2 = new Particle(Pion2, zero, rcol
369 pion3 = new Particle(Pion3, zero, rcol
370 } else if (rdm < 2.0 / 3.0) {
371 pion1 = new Particle(Pion1, zero, rcol
372 pion2 = new Particle(Pion3, zero, rcol
373 pion3 = new Particle(Pion2, zero, rcol
374 } else {
375 pion1 = new Particle(Pion2, zero, rcol
376 pion2 = new Particle(Pion1, zero, rcol
377 pion3 = new Particle(Pion3, zero, rcol
378 }
379
380 list.push_back(pion1);
381 list.push_back(pion2);
382 list.push_back(pion3);
383
384 PhaseSpaceGenerator::generate(sqrtS, list)
385
386 fs->addModifiedParticle(nucleon);
387 fs->addModifiedParticle(antinucleon);
388 fs->addCreatedParticle(pion1);
389 fs->addCreatedParticle(pion2);
390 fs->addCreatedParticle(pion3);
391
392 }
393 }
394