| 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 "G4INCLNNbarToNNbar2piChannel.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 NNbarToNNbar2piChannel::NNbarToNNbar2piChannel(Particle *p1, Particle *p2)
50 : particle1(p1), particle2(p2)
51 {}
52
53 NNbarToNNbar2piChannel::~NNbarToNNbar2piChannel(){}
54
55 void NNbarToNNbar2piChannel::fillFinalState(FinalState *fs) {
56
57 //brief ppbar
58 // p pbar -> p pbar pi+ pi- (BFMM 167)
59 // p pbar -> p nbar pi- pi0 (same as BFMM 490)
60 // p pbar -> n pbar pi+ pi0 (same as BFMM 490)
61 // p pbar -> n nbar pi+ pi- (BFMM 198)
62 //
63 //brief npbar
64 // n pbar -> p pbar pi- pi0 (BFMM 490)
65 // n pbar -> p nbar pi- pi- (BFMM 492)
66 // n pbar -> n nbar pi- pi0 (same as BFMM 490)
67 // n pbar -> n pbar pi+ pi- (BFMM 494)
68 //
69 //brief nnbar
70 // n nbar -> n nbar pi+ pi- (same as BFMM 167)
71 // n nbar -> p nbar pi- pi0 (same as BFMM 490)
72 // n nbar -> n pbar pi+ pi0 (same as BFMM 490)
73 // n nbar -> p pbar pi+ pi- (same as BFMM 198)
74 //
75 //brief pnbar
76 // p nbar -> p pbar pi+ pi0 (same as BFMM 490)
77 // p nbar -> n pbar pi+ pi+ (same as BFMM 492)
78 // p nbar -> n nbar pi+ pi0 (same as BFMM 490)
79 // p nbar -> p nbar pi+ pi- (same as BFMM 494)
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); //GeV
94 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(nucleon, antinucleon);
95 const G4double rdm = Random::shoot();
96
97 const std::vector<G4double> BFMM167 = {-6.885, 0.476, 1.206, 13.857, -5.728, 1.220};
98 //const G4double Eth_PPbar_PPbar_pip_pim = 1.220;
99 const std::vector<G4double> BFMM198 = {1.857, -21.213, -3.448, 0.827, -0.390, 1.231};
100 //const G4double Eth_PPbar_NNbar_pip_pim = 1.231;
101 const std::vector<G4double> BFMM490 = {-3.594, 0.811, 0.306, 5.108, -1.625, 1.201};
102 //const G4double Eth_PNbar_PPbar_pim_pi0 = 1.201;
103 const std::vector<G4double> BFMM492 = {-5.443, 7.254, -2.936, 8.441, -2.588, 1.221};
104 //const G4double Eth_PNbar_NPbar_pim_pim = 1.221;
105 const std::vector<G4double> BFMM494 = {21.688, -38.709, -2.062, -17.783, 3.895, 1.221};
106 //const G4double Eth_NPbar_NPbar_pip_pim = 1.221;
107
108 // pnbar total is same as for npbar
109 // ppbar total is same as for nnbar
110 const G4double totalppbar = KinematicsUtils::compute_xs(BFMM167, plab) +KinematicsUtils::compute_xs(std::move(BFMM198), plab) +2*KinematicsUtils::compute_xs(BFMM490, plab);
111 const G4double totalpnbar = KinematicsUtils::compute_xs(BFMM492, plab) +KinematicsUtils::compute_xs(std::move(BFMM494), plab) +2*KinematicsUtils::compute_xs(BFMM490, plab);
112 //totalnnbar == totalppbar;
113 //totalpnbar == totalnpbar;
114 ParticleType Pion1;
115 ParticleType Pion2;
116
117 //setting types of new particles
118 if(nucleon->getType()==Proton){
119 if(antinucleon->getType()==antiProton){ // ppbar case
120 if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)){ // ppbarpi-pi+ case
121 Pion1 = PiMinus;
122 Pion2 = PiPlus;
123 if(rdm<0.5){
124 nucleon->setType(Proton);
125 antinucleon->setType(antiProton);
126 }
127 else{
128 nucleon->setType(antiProton);
129 antinucleon->setType(Proton);
130 }
131 }
132 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)+KinematicsUtils::compute_xs(BFMM490, plab)){ //pnbarpi-pi0 case
133 Pion1 = PiMinus;
134 Pion2 = PiZero;
135 if(rdm<0.5){
136 nucleon->setType(Proton);
137 antinucleon->setType(antiNeutron);
138 }
139 else{
140 nucleon->setType(antiNeutron);
141 antinucleon->setType(Proton);
142 }
143 }
144 else if(rdm*totalppbar < KinematicsUtils::compute_xs(std::move(BFMM167), plab)+2*KinematicsUtils::compute_xs(std::move(BFMM490), plab)){ //npbarpi+pi0 case
145 Pion1 = PiPlus;
146 Pion2 = PiZero;
147 if(rdm<0.5){
148 nucleon->setType(Neutron);
149 antinucleon->setType(antiProton);
150 }
151 else{
152 nucleon->setType(antiProton);
153 antinucleon->setType(Neutron);
154 }
155 }
156 else{ // n nbar pi+ pi- case case
157 Pion1 = PiMinus;
158 Pion2 = PiPlus;
159 if(rdm<0.5){
160 nucleon->setType(Neutron);
161 antinucleon->setType(antiNeutron);
162 }
163 else{
164 nucleon->setType(antiNeutron);
165 antinucleon->setType(Neutron);
166 }
167 }
168 }
169 else{ //antiNeutron (pnbar case)
170 if(rdm*totalpnbar < KinematicsUtils::compute_xs(BFMM490, plab)){ // p pbar pi+ pi0 case
171 Pion1 = PiZero;
172 Pion2 = PiPlus;
173 if(rdm<0.5){
174 nucleon->setType(Proton);
175 antinucleon->setType(antiProton);
176 }
177 else{
178 nucleon->setType(antiProton);
179 antinucleon->setType(Proton);
180 }
181 }
182 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM490, plab)+KinematicsUtils::compute_xs(BFMM492, plab)){ // n pbar pi+ pi+ case
183 Pion1 = PiPlus;
184 Pion2 = PiPlus;
185 if(rdm<0.5){
186 nucleon->setType(Neutron);
187 antinucleon->setType(antiProton);
188 }
189 else{
190 nucleon->setType(antiProton);
191 antinucleon->setType(Neutron);
192 }
193 }
194 else if(rdm*totalppbar < 2*KinematicsUtils::compute_xs(std::move(BFMM490), plab)+KinematicsUtils::compute_xs(std::move(BFMM492), plab)){ // n nbar pi+ pi0 case
195 Pion1 = PiZero;
196 Pion2 = PiPlus;
197 if(rdm<0.5){
198 nucleon->setType(Neutron);
199 antinucleon->setType(antiNeutron);
200 }
201 else{
202 nucleon->setType(antiNeutron);
203 antinucleon->setType(Neutron);
204 }
205 }
206 else{ // p nbar pi+ pi- case
207 Pion1 = PiMinus;
208 Pion2 = PiPlus;
209 if(rdm<0.5){
210 nucleon->setType(Proton);
211 antinucleon->setType(antiNeutron);
212 }
213 else{
214 nucleon->setType(antiNeutron);
215 antinucleon->setType(Proton);
216 }
217 }
218 }
219 }
220 else{ // neutron
221 if(antinucleon->getType()==antiProton){ //npbar case
222 if(rdm*totalpnbar < KinematicsUtils::compute_xs(BFMM490, plab)){ // p pbar pi- pi0 case
223 Pion1 = PiZero;
224 Pion2 = PiMinus;
225 if(rdm<0.5){
226 nucleon->setType(Proton);
227 antinucleon->setType(antiProton);
228 }
229 else{
230 nucleon->setType(antiProton);
231 antinucleon->setType(Proton);
232 }
233 }
234 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM490, plab)+KinematicsUtils::compute_xs(BFMM492, plab)){ // p nbar pi- pi- case
235 Pion1 = PiMinus;
236 Pion2 = PiMinus;
237 if(rdm<0.5){
238 nucleon->setType(Proton);
239 antinucleon->setType(antiNeutron);
240 }
241 else{
242 nucleon->setType(antiNeutron);
243 antinucleon->setType(Proton);
244 }
245 }
246 else if(rdm*totalppbar < 2*KinematicsUtils::compute_xs(std::move(BFMM490), plab)+KinematicsUtils::compute_xs(std::move(BFMM492), plab)){ // n nbar pi- pi0 case
247 Pion1 = PiZero;
248 Pion2 = PiMinus;
249 if(rdm<0.5){
250 nucleon->setType(Neutron);
251 antinucleon->setType(antiNeutron);
252 }
253 else{
254 nucleon->setType(antiNeutron);
255 antinucleon->setType(Neutron);
256 }
257 }
258 else{ // n pbar pi+ pi- case
259 Pion1 = PiMinus;
260 Pion2 = PiPlus;
261 if(rdm<0.5){
262 nucleon->setType(Neutron);
263 antinucleon->setType(antiProton);
264 }
265 else{
266 nucleon->setType(antiProton);
267 antinucleon->setType(Neutron);
268 }
269 }
270 }
271 else{ //antiNeutron (nnbar case)
272 if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)){ // nnbarpi-pi+ case
273 Pion1 = PiMinus;
274 Pion2 = PiPlus;
275 if(rdm<0.5){
276 nucleon->setType(Neutron);
277 antinucleon->setType(antiNeutron);
278 }
279 else{
280 nucleon->setType(antiNeutron);
281 antinucleon->setType(Neutron);
282 }
283 }
284 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)+KinematicsUtils::compute_xs(BFMM490, plab)){ //pnbarpi-pi0 case
285 Pion1 = PiMinus;
286 Pion2 = PiZero;
287 if(rdm<0.5){
288 nucleon->setType(Proton);
289 antinucleon->setType(antiNeutron);
290 }
291 else{
292 nucleon->setType(antiNeutron);
293 antinucleon->setType(Proton);
294 }
295 }
296 else if(rdm*totalppbar < KinematicsUtils::compute_xs(std::move(BFMM167), plab)+2*KinematicsUtils::compute_xs(std::move(BFMM490), plab)){ //npbarpi+pi0 case
297 Pion1 = PiPlus;
298 Pion2 = PiZero;
299 if(rdm<0.5){
300 nucleon->setType(Neutron);
301 antinucleon->setType(antiProton);
302 }
303 else{
304 nucleon->setType(antiProton);
305 antinucleon->setType(Neutron);
306 }
307 }
308 else{ // p pbar pi+ pi- case
309 Pion1 = PiMinus;
310 Pion2 = PiPlus;
311 if(rdm<0.5){
312 nucleon->setType(Proton);
313 antinucleon->setType(antiProton);
314 }
315 else{
316 nucleon->setType(antiProton);
317 antinucleon->setType(Proton);
318 }
319 }
320 }
321 }
322
323 ParticleList list;
324 list.push_back(nucleon);
325 list.push_back(antinucleon);
326 const ThreeVector &rcol = nucleon->getPosition();
327 const ThreeVector zero;
328
329 Particle *pion2 = new Particle(Pion1,zero,rcol);
330 Particle *pion1 = new Particle(Pion2,zero,rcol);
331 if(rdm < 0.5){
332 pion1->setType(Pion1);
333 pion2->setType(Pion2);
334 }
335
336 list.push_back(pion1);
337 list.push_back(pion2);
338
339 PhaseSpaceGenerator::generate(sqrtS, list);
340
341 fs->addModifiedParticle(nucleon);
342 fs->addModifiedParticle(antinucleon);
343 fs->addCreatedParticle(pion1);
344 fs->addCreatedParticle(pion2);
345
346 }
347 }
348