| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 #include "G4INCLNNbarToNNbar2piChannel.hh" 38 #include "G4INCLNNbarToNNbar2piChannel.hh"
39 #include "G4INCLKinematicsUtils.hh" 39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh" 40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh" 41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh" 42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh" 43 #include "G4INCLLogger.hh"
44 #include <algorithm> 44 #include <algorithm>
45 #include "G4INCLPhaseSpaceGenerator.hh" 45 #include "G4INCLPhaseSpaceGenerator.hh"
46 46
47 namespace G4INCL { 47 namespace G4INCL {
48 48
49 NNbarToNNbar2piChannel::NNbarToNNbar2piChann 49 NNbarToNNbar2piChannel::NNbarToNNbar2piChannel(Particle *p1, Particle *p2)
50 : particle1(p1), particle2(p2) 50 : particle1(p1), particle2(p2)
51 {} 51 {}
52 52
53 NNbarToNNbar2piChannel::~NNbarToNNbar2piChan 53 NNbarToNNbar2piChannel::~NNbarToNNbar2piChannel(){}
54 54
55 void NNbarToNNbar2piChannel::fillFinalState( 55 void NNbarToNNbar2piChannel::fillFinalState(FinalState *fs) {
56 56
57 //brief ppbar 57 //brief ppbar
58 // p pbar -> p pbar pi+ pi- (BFMM 167) 58 // p pbar -> p pbar pi+ pi- (BFMM 167)
59 // p pbar -> p nbar pi- pi0 (same as B 59 // p pbar -> p nbar pi- pi0 (same as BFMM 490)
60 // p pbar -> n pbar pi+ pi0 (same as B 60 // p pbar -> n pbar pi+ pi0 (same as BFMM 490)
61 // p pbar -> n nbar pi+ pi- (BFMM 198) 61 // p pbar -> n nbar pi+ pi- (BFMM 198)
62 // 62 //
63 //brief npbar 63 //brief npbar
64 // n pbar -> p pbar pi- pi0 (BFMM 490) 64 // n pbar -> p pbar pi- pi0 (BFMM 490)
65 // n pbar -> p nbar pi- pi- (BFMM 492) 65 // n pbar -> p nbar pi- pi- (BFMM 492)
66 // n pbar -> n nbar pi- pi0 (same as B 66 // n pbar -> n nbar pi- pi0 (same as BFMM 490)
67 // n pbar -> n pbar pi+ pi- (BFMM 494) 67 // n pbar -> n pbar pi+ pi- (BFMM 494)
68 // 68 //
69 //brief nnbar 69 //brief nnbar
70 // n nbar -> n nbar pi+ pi- (same as B 70 // n nbar -> n nbar pi+ pi- (same as BFMM 167)
71 // n nbar -> p nbar pi- pi0 (same as B 71 // n nbar -> p nbar pi- pi0 (same as BFMM 490)
72 // n nbar -> n pbar pi+ pi0 (same as B 72 // n nbar -> n pbar pi+ pi0 (same as BFMM 490)
73 // n nbar -> p pbar pi+ pi- (same as B 73 // n nbar -> p pbar pi+ pi- (same as BFMM 198)
74 // 74 //
75 //brief pnbar 75 //brief pnbar
76 // p nbar -> p pbar pi+ pi0 (same as B 76 // p nbar -> p pbar pi+ pi0 (same as BFMM 490)
77 // p nbar -> n pbar pi+ pi+ (same as B 77 // p nbar -> n pbar pi+ pi+ (same as BFMM 492)
78 // p nbar -> n nbar pi+ pi0 (same as B 78 // p nbar -> n nbar pi+ pi0 (same as BFMM 490)
79 // p nbar -> p nbar pi+ pi- (same as B 79 // p nbar -> p nbar pi+ pi- (same as BFMM 494)
80 80
81 Particle *nucleon; 81 Particle *nucleon;
82 Particle *antinucleon; 82 Particle *antinucleon;
83 83
84 if(particle1->isNucleon()){ 84 if(particle1->isNucleon()){
85 nucleon = particle1; 85 nucleon = particle1;
86 antinucleon = particle2; 86 antinucleon = particle2;
87 } 87 }
88 else{ 88 else{
89 nucleon = particle2; 89 nucleon = particle2;
90 antinucleon = particle1; 90 antinucleon = particle1;
91 } 91 }
92 92
93 const G4double plab = 0.001*KinematicsUtil 93 const G4double plab = 0.001*KinematicsUtils::momentumInLab(particle1, particle2); //GeV
94 const G4double sqrtS = KinematicsUtils::to 94 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(nucleon, antinucleon);
95 const G4double rdm = Random::shoot(); 95 const G4double rdm = Random::shoot();
96 96
97 const std::vector<G4double> BFMM167 = 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_p 98 //const G4double Eth_PPbar_PPbar_pip_pim = 1.220;
99 const std::vector<G4double> BFMM198 = 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_p 100 //const G4double Eth_PPbar_NNbar_pip_pim = 1.231;
101 const std::vector<G4double> BFMM490 = 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_p 102 //const G4double Eth_PNbar_PPbar_pim_pi0 = 1.201;
103 const std::vector<G4double> BFMM492 = 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_p 104 //const G4double Eth_PNbar_NPbar_pim_pim = 1.221;
105 const std::vector<G4double> BFMM494 = 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_p 106 //const G4double Eth_NPbar_NPbar_pip_pim = 1.221;
107 107
108 // pnbar total is same as for npbar 108 // pnbar total is same as for npbar
109 // ppbar total is same as for nnbar 109 // ppbar total is same as for nnbar
110 const G4double totalppbar = KinematicsUtil << 110 const G4double totalppbar = KinematicsUtils::compute_xs(BFMM167, plab) +KinematicsUtils::compute_xs(BFMM198, plab) +2*KinematicsUtils::compute_xs(BFMM490, plab);
111 const G4double totalpnbar = KinematicsUtil << 111 const G4double totalpnbar = KinematicsUtils::compute_xs(BFMM492, plab) +KinematicsUtils::compute_xs(BFMM494, plab) +2*KinematicsUtils::compute_xs(BFMM490, plab);
112 //totalnnbar == totalppbar; 112 //totalnnbar == totalppbar;
113 //totalpnbar == totalnpbar; 113 //totalpnbar == totalnpbar;
114 ParticleType Pion1; 114 ParticleType Pion1;
115 ParticleType Pion2; 115 ParticleType Pion2;
116 116
117 //setting types of new particles 117 //setting types of new particles
118 if(nucleon->getType()==Proton){ 118 if(nucleon->getType()==Proton){
119 if(antinucleon->getType()==antiProton){ 119 if(antinucleon->getType()==antiProton){ // ppbar case
120 if(rdm*totalppbar < KinematicsUtils::c 120 if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)){ // ppbarpi-pi+ case
121 Pion1 = PiMinus; 121 Pion1 = PiMinus;
122 Pion2 = PiPlus; 122 Pion2 = PiPlus;
123 if(rdm<0.5){ 123 if(rdm<0.5){
124 nucleon->setType(Proton); 124 nucleon->setType(Proton);
125 antinucleon->setType(antiProton); 125 antinucleon->setType(antiProton);
126 } 126 }
127 else{ 127 else{
128 nucleon->setType(antiProton); 128 nucleon->setType(antiProton);
129 antinucleon->setType(Proton); 129 antinucleon->setType(Proton);
130 } 130 }
131 } 131 }
132 else if(rdm*totalppbar < KinematicsUti 132 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)+KinematicsUtils::compute_xs(BFMM490, plab)){ //pnbarpi-pi0 case
133 Pion1 = PiMinus; 133 Pion1 = PiMinus;
134 Pion2 = PiZero; 134 Pion2 = PiZero;
135 if(rdm<0.5){ 135 if(rdm<0.5){
136 nucleon->setType(Proton); 136 nucleon->setType(Proton);
137 antinucleon->setType(antiNeutron); 137 antinucleon->setType(antiNeutron);
138 } 138 }
139 else{ 139 else{
140 nucleon->setType(antiNeutron); 140 nucleon->setType(antiNeutron);
141 antinucleon->setType(Proton); 141 antinucleon->setType(Proton);
142 } 142 }
143 } 143 }
144 else if(rdm*totalppbar < KinematicsUti << 144 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)+2*KinematicsUtils::compute_xs(BFMM490, plab)){ //npbarpi+pi0 case
145 Pion1 = PiPlus; 145 Pion1 = PiPlus;
146 Pion2 = PiZero; 146 Pion2 = PiZero;
147 if(rdm<0.5){ 147 if(rdm<0.5){
148 nucleon->setType(Neutron); 148 nucleon->setType(Neutron);
149 antinucleon->setType(antiProton); 149 antinucleon->setType(antiProton);
150 } 150 }
151 else{ 151 else{
152 nucleon->setType(antiProton); 152 nucleon->setType(antiProton);
153 antinucleon->setType(Neutron); 153 antinucleon->setType(Neutron);
154 } 154 }
155 } 155 }
156 else{ // n nbar pi+ pi- case case 156 else{ // n nbar pi+ pi- case case
157 Pion1 = PiMinus; 157 Pion1 = PiMinus;
158 Pion2 = PiPlus; 158 Pion2 = PiPlus;
159 if(rdm<0.5){ 159 if(rdm<0.5){
160 nucleon->setType(Neutron); 160 nucleon->setType(Neutron);
161 antinucleon->setType(antiNeutron); 161 antinucleon->setType(antiNeutron);
162 } 162 }
163 else{ 163 else{
164 nucleon->setType(antiNeutron); 164 nucleon->setType(antiNeutron);
165 antinucleon->setType(Neutron); 165 antinucleon->setType(Neutron);
166 } 166 }
167 } 167 }
168 } 168 }
169 else{ //antiNeutron (pnbar case) 169 else{ //antiNeutron (pnbar case)
170 if(rdm*totalpnbar < KinematicsUtils::c 170 if(rdm*totalpnbar < KinematicsUtils::compute_xs(BFMM490, plab)){ // p pbar pi+ pi0 case
171 Pion1 = PiZero; 171 Pion1 = PiZero;
172 Pion2 = PiPlus; 172 Pion2 = PiPlus;
173 if(rdm<0.5){ 173 if(rdm<0.5){
174 nucleon->setType(Proton); 174 nucleon->setType(Proton);
175 antinucleon->setType(antiProton); 175 antinucleon->setType(antiProton);
176 } 176 }
177 else{ 177 else{
178 nucleon->setType(antiProton); 178 nucleon->setType(antiProton);
179 antinucleon->setType(Proton); 179 antinucleon->setType(Proton);
180 } 180 }
181 } 181 }
182 else if(rdm*totalppbar < KinematicsUti 182 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM490, plab)+KinematicsUtils::compute_xs(BFMM492, plab)){ // n pbar pi+ pi+ case
183 Pion1 = PiPlus; 183 Pion1 = PiPlus;
184 Pion2 = PiPlus; 184 Pion2 = PiPlus;
185 if(rdm<0.5){ 185 if(rdm<0.5){
186 nucleon->setType(Neutron); 186 nucleon->setType(Neutron);
187 antinucleon->setType(antiProton); 187 antinucleon->setType(antiProton);
188 } 188 }
189 else{ 189 else{
190 nucleon->setType(antiProton); 190 nucleon->setType(antiProton);
191 antinucleon->setType(Neutron); 191 antinucleon->setType(Neutron);
192 } 192 }
193 } 193 }
194 else if(rdm*totalppbar < 2*KinematicsU << 194 else if(rdm*totalppbar < 2*KinematicsUtils::compute_xs(BFMM490, plab)+KinematicsUtils::compute_xs(BFMM492, plab)){ // n nbar pi+ pi0 case
195 Pion1 = PiZero; 195 Pion1 = PiZero;
196 Pion2 = PiPlus; 196 Pion2 = PiPlus;
197 if(rdm<0.5){ 197 if(rdm<0.5){
198 nucleon->setType(Neutron); 198 nucleon->setType(Neutron);
199 antinucleon->setType(antiNeutron); 199 antinucleon->setType(antiNeutron);
200 } 200 }
201 else{ 201 else{
202 nucleon->setType(antiNeutron); 202 nucleon->setType(antiNeutron);
203 antinucleon->setType(Neutron); 203 antinucleon->setType(Neutron);
204 } 204 }
205 } 205 }
206 else{ // p nbar pi+ pi- case 206 else{ // p nbar pi+ pi- case
207 Pion1 = PiMinus; 207 Pion1 = PiMinus;
208 Pion2 = PiPlus; 208 Pion2 = PiPlus;
209 if(rdm<0.5){ 209 if(rdm<0.5){
210 nucleon->setType(Proton); 210 nucleon->setType(Proton);
211 antinucleon->setType(antiNeutron); 211 antinucleon->setType(antiNeutron);
212 } 212 }
213 else{ 213 else{
214 nucleon->setType(antiNeutron); 214 nucleon->setType(antiNeutron);
215 antinucleon->setType(Proton); 215 antinucleon->setType(Proton);
216 } 216 }
217 } 217 }
218 } 218 }
219 } 219 }
220 else{ // neutron 220 else{ // neutron
221 if(antinucleon->getType()==antiProton){ 221 if(antinucleon->getType()==antiProton){ //npbar case
222 if(rdm*totalpnbar < KinematicsUtils::c 222 if(rdm*totalpnbar < KinematicsUtils::compute_xs(BFMM490, plab)){ // p pbar pi- pi0 case
223 Pion1 = PiZero; 223 Pion1 = PiZero;
224 Pion2 = PiMinus; 224 Pion2 = PiMinus;
225 if(rdm<0.5){ 225 if(rdm<0.5){
226 nucleon->setType(Proton); 226 nucleon->setType(Proton);
227 antinucleon->setType(antiProton); 227 antinucleon->setType(antiProton);
228 } 228 }
229 else{ 229 else{
230 nucleon->setType(antiProton); 230 nucleon->setType(antiProton);
231 antinucleon->setType(Proton); 231 antinucleon->setType(Proton);
232 } 232 }
233 } 233 }
234 else if(rdm*totalppbar < KinematicsUti 234 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM490, plab)+KinematicsUtils::compute_xs(BFMM492, plab)){ // p nbar pi- pi- case
235 Pion1 = PiMinus; 235 Pion1 = PiMinus;
236 Pion2 = PiMinus; 236 Pion2 = PiMinus;
237 if(rdm<0.5){ 237 if(rdm<0.5){
238 nucleon->setType(Proton); 238 nucleon->setType(Proton);
239 antinucleon->setType(antiNeutron); 239 antinucleon->setType(antiNeutron);
240 } 240 }
241 else{ 241 else{
242 nucleon->setType(antiNeutron); 242 nucleon->setType(antiNeutron);
243 antinucleon->setType(Proton); 243 antinucleon->setType(Proton);
244 } 244 }
245 } 245 }
246 else if(rdm*totalppbar < 2*KinematicsU << 246 else if(rdm*totalppbar < 2*KinematicsUtils::compute_xs(BFMM490, plab)+KinematicsUtils::compute_xs(BFMM492, plab)){ // n nbar pi- pi0 case
247 Pion1 = PiZero; 247 Pion1 = PiZero;
248 Pion2 = PiMinus; 248 Pion2 = PiMinus;
249 if(rdm<0.5){ 249 if(rdm<0.5){
250 nucleon->setType(Neutron); 250 nucleon->setType(Neutron);
251 antinucleon->setType(antiNeutron); 251 antinucleon->setType(antiNeutron);
252 } 252 }
253 else{ 253 else{
254 nucleon->setType(antiNeutron); 254 nucleon->setType(antiNeutron);
255 antinucleon->setType(Neutron); 255 antinucleon->setType(Neutron);
256 } 256 }
257 } 257 }
258 else{ // n pbar pi+ pi- case 258 else{ // n pbar pi+ pi- case
259 Pion1 = PiMinus; 259 Pion1 = PiMinus;
260 Pion2 = PiPlus; 260 Pion2 = PiPlus;
261 if(rdm<0.5){ 261 if(rdm<0.5){
262 nucleon->setType(Neutron); 262 nucleon->setType(Neutron);
263 antinucleon->setType(antiProton); 263 antinucleon->setType(antiProton);
264 } 264 }
265 else{ 265 else{
266 nucleon->setType(antiProton); 266 nucleon->setType(antiProton);
267 antinucleon->setType(Neutron); 267 antinucleon->setType(Neutron);
268 } 268 }
269 } 269 }
270 } 270 }
271 else{ //antiNeutron (nnbar case) 271 else{ //antiNeutron (nnbar case)
272 if(rdm*totalppbar < KinematicsUtils::c 272 if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)){ // nnbarpi-pi+ case
273 Pion1 = PiMinus; 273 Pion1 = PiMinus;
274 Pion2 = PiPlus; 274 Pion2 = PiPlus;
275 if(rdm<0.5){ 275 if(rdm<0.5){
276 nucleon->setType(Neutron); 276 nucleon->setType(Neutron);
277 antinucleon->setType(antiNeutron); 277 antinucleon->setType(antiNeutron);
278 } 278 }
279 else{ 279 else{
280 nucleon->setType(antiNeutron); 280 nucleon->setType(antiNeutron);
281 antinucleon->setType(Neutron); 281 antinucleon->setType(Neutron);
282 } 282 }
283 } 283 }
284 else if(rdm*totalppbar < KinematicsUti 284 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)+KinematicsUtils::compute_xs(BFMM490, plab)){ //pnbarpi-pi0 case
285 Pion1 = PiMinus; 285 Pion1 = PiMinus;
286 Pion2 = PiZero; 286 Pion2 = PiZero;
287 if(rdm<0.5){ 287 if(rdm<0.5){
288 nucleon->setType(Proton); 288 nucleon->setType(Proton);
289 antinucleon->setType(antiNeutron); 289 antinucleon->setType(antiNeutron);
290 } 290 }
291 else{ 291 else{
292 nucleon->setType(antiNeutron); 292 nucleon->setType(antiNeutron);
293 antinucleon->setType(Proton); 293 antinucleon->setType(Proton);
294 } 294 }
295 } 295 }
296 else if(rdm*totalppbar < KinematicsUti << 296 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM167, plab)+2*KinematicsUtils::compute_xs(BFMM490, plab)){ //npbarpi+pi0 case
297 Pion1 = PiPlus; 297 Pion1 = PiPlus;
298 Pion2 = PiZero; 298 Pion2 = PiZero;
299 if(rdm<0.5){ 299 if(rdm<0.5){
300 nucleon->setType(Neutron); 300 nucleon->setType(Neutron);
301 antinucleon->setType(antiProton); 301 antinucleon->setType(antiProton);
302 } 302 }
303 else{ 303 else{
304 nucleon->setType(antiProton); 304 nucleon->setType(antiProton);
305 antinucleon->setType(Neutron); 305 antinucleon->setType(Neutron);
306 } 306 }
307 } 307 }
308 else{ // p pbar pi+ pi- case 308 else{ // p pbar pi+ pi- case
309 Pion1 = PiMinus; 309 Pion1 = PiMinus;
310 Pion2 = PiPlus; 310 Pion2 = PiPlus;
311 if(rdm<0.5){ 311 if(rdm<0.5){
312 nucleon->setType(Proton); 312 nucleon->setType(Proton);
313 antinucleon->setType(antiProton); 313 antinucleon->setType(antiProton);
314 } 314 }
315 else{ 315 else{
316 nucleon->setType(antiProton); 316 nucleon->setType(antiProton);
317 antinucleon->setType(Proton); 317 antinucleon->setType(Proton);
318 } 318 }
319 } 319 }
320 } 320 }
321 } 321 }
322 322
323 ParticleList list; 323 ParticleList list;
324 list.push_back(nucleon); 324 list.push_back(nucleon);
325 list.push_back(antinucleon); 325 list.push_back(antinucleon);
326 const ThreeVector &rcol = nucleon->getPosi 326 const ThreeVector &rcol = nucleon->getPosition();
327 const ThreeVector zero; 327 const ThreeVector zero;
328 328
329 Particle *pion2 = new Particle(Pion1,zero, 329 Particle *pion2 = new Particle(Pion1,zero,rcol);
330 Particle *pion1 = new Particle(Pion2,zero, 330 Particle *pion1 = new Particle(Pion2,zero,rcol);
331 if(rdm < 0.5){ 331 if(rdm < 0.5){
332 pion1->setType(Pion1); 332 pion1->setType(Pion1);
333 pion2->setType(Pion2); 333 pion2->setType(Pion2);
334 } 334 }
335 335
336 list.push_back(pion1); 336 list.push_back(pion1);
337 list.push_back(pion2); 337 list.push_back(pion2);
338 338
339 PhaseSpaceGenerator::generate(sqrtS, list) 339 PhaseSpaceGenerator::generate(sqrtS, list);
340 340
341 fs->addModifiedParticle(nucleon); 341 fs->addModifiedParticle(nucleon);
342 fs->addModifiedParticle(antinucleon); 342 fs->addModifiedParticle(antinucleon);
343 fs->addCreatedParticle(pion1); 343 fs->addCreatedParticle(pion1);
344 fs->addCreatedParticle(pion2); 344 fs->addCreatedParticle(pion2);
345 345
346 } 346 }
347 } 347 }
348 348