| Geant4 Cross Reference |
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
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::NNbarToNNbar3piChannel(Particle *p1, Particle *p2)
50 : particle1(p1), particle2(p2)
51 {}
52
53 NNbarToNNbar3piChannel::~NNbarToNNbar3piChannel(){}
54
55 void NNbarToNNbar3piChannel::fillFinalState(FinalState *fs) {
56
57 //brief ppbar
58 // p pbar -> p pbar pi+ pi- pi0 (BFMM 161)
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 197)
62 //
63 //brief npbar
64 // n pbar -> p pbar 2pi- pi+ (same as BFMM 169)
65 // n pbar -> p nbar 2pi- pi0 (same as BFMM 197)
66 // n pbar -> n nbar 2pi- pi+ (same as BFMM 169)
67 // n pbar -> n pbar pi+ pi- pi0 (same as BFMM 161)
68 //
69 //brief nnbar
70 // n nbar -> n nbar pi+ pi- pi0 (same as BFMM 161)
71 // n nbar -> p nbar 2pi- pi+ (same as BFMM 169)
72 // n nbar -> n pbar 2pi+ pi- (same as BFMM 201)
73 // n nbar -> p pbar pi+ pi- pi0 (same as BFMM 197)
74 //
75 //brief pnbar
76 // p nbar -> p pbar 2pi+ pi- (same as BFMM 169)
77 // p nbar -> n pbar 2pi+ pi0 (same as BFMM 197)
78 // p nbar -> n nbar 2pi+ pi- (same as BFMM 169)
79 // p nbar -> p nbar pi+ pi- pi0 (same as BFMM 161)
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*KinematicsUtils::momentumInLab(particle1, particle2);
94 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(nucleon, antinucleon);
95 const G4double rdm = Random::shoot();
96
97 const std::vector<G4double> BFMM161 = {-6.434, 1.351, -5.185, 7.754, -1.692, 1.604};
98 //const G4double Eth_PPbar_PPbar_pip_pim_pi0 = 1.604;
99 const std::vector<G4double> BFMM169 = {3.696, -5.356, -0.053, 1.941, -0.432, 1.624};
100 //const G4double Eth_PPbar_PNbar_2pim_pip = 1.624;
101 const std::vector<G4double> BFMM201 = {-1.070, -0.636, -0.009, 2.335, -0.499, 1.624};
102 //const G4double Eth_PPbar_NPbar_2pip_pim = 1.624;
103 const std::vector<G4double> BFMM197 = {1.857, -21.213, -3.448, 0.827, -0.390, 1.616};
104 //const G4double Eth_PPbar_NNbar_pip_pim_pi0 = 1.616;
105
106 // pnbar total is same as for npbar
107 // ppbar total is same as for nnbar
108 const G4double totalppbar = KinematicsUtils::compute_xs(BFMM161, plab)
109 +KinematicsUtils::compute_xs(BFMM169, plab)
110 +KinematicsUtils::compute_xs(BFMM201, plab)
111 +KinematicsUtils::compute_xs(BFMM197, plab);
112 const G4double totalpnbar = KinematicsUtils::compute_xs(BFMM161, plab)
113 +KinematicsUtils::compute_xs(BFMM197, plab)
114 +2*KinematicsUtils::compute_xs(BFMM169, plab);
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){ // ppbar case
125 if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM161, plab)){ // p pbar pi+ pi- pi0 case
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 < KinematicsUtils::compute_xs(BFMM161, plab)
139 +KinematicsUtils::compute_xs(BFMM169, plab)){ //p nbar 2pi- pi+ case
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 < KinematicsUtils::compute_xs(std::move(BFMM161), plab)
153 +KinematicsUtils::compute_xs(std::move(BFMM169), plab)
154 +KinematicsUtils::compute_xs(std::move(BFMM201), plab)){ //n pbar 2pi+ pi- case
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::compute_xs(BFMM169, plab)){ // p pbar 2pi+ pi- case
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 < KinematicsUtils::compute_xs(BFMM169, plab)
196 +KinematicsUtils::compute_xs(BFMM197, plab)){ // n pbar 2pi+ pi0 case
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*KinematicsUtils::compute_xs(std::move(BFMM169), plab)
210 +KinematicsUtils::compute_xs(std::move(BFMM197), plab)){ // n nbar 2pi+ pi- case
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){ //npbar case
240 if(rdm*totalpnbar < KinematicsUtils::compute_xs(BFMM169, plab)){ // p pbar 2pi- pi+ case
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 < KinematicsUtils::compute_xs(BFMM169, plab)
254 +KinematicsUtils::compute_xs(BFMM197, plab)){ // p nbar 2pi- pi0 case
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*KinematicsUtils::compute_xs(std::move(BFMM169), plab)
268 +KinematicsUtils::compute_xs(std::move(BFMM197), plab)){ // n nbar 2pi- pi+ case
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::compute_xs(BFMM161, plab)){ // n nbar pi+ pi- pi0 case
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 < KinematicsUtils::compute_xs(BFMM161, plab)
310 +KinematicsUtils::compute_xs(BFMM169, plab)){ //p nbar 2pi- pi+ case
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 < KinematicsUtils::compute_xs(std::move(BFMM161), plab)
324 +KinematicsUtils::compute_xs(std::move(BFMM169), plab)
325 +KinematicsUtils::compute_xs(std::move(BFMM201), plab)){ //n pbar 2pi+ pi- case
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->getPosition();
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 on the random number
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