Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCLNNbarToAnnihilationChannel.cc
Go to the documentation of this file.
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#include "G4EnvironmentUtils.hh"
38#include "G4INCLNNbarToAnnihilationChannel.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
47namespace G4INCL {
48
49 NNbarToAnnihilationChannel::NNbarToAnnihilationChannel(Nucleus *n, Particle *p1, Particle *p2)
50 :theNucleus(n), particle1(p1), particle2(p2)
51 {}
52
53 NNbarToAnnihilationChannel::~NNbarToAnnihilationChannel(){}
54
55//fill probabilities and particle types from datafile and return probability sum for normalization
56G4double NNbarToAnnihilationChannel::read_file(std::string filename, std::vector<G4double>& probabilities, std::vector<std::vector<std::string>>& particle_types) {
57 std::ifstream file(filename);
58 G4double sum_probs = 0.0;
59 if (file.is_open()) {
60 std::string line;
61 while (getline(file, line)) {
62 std::istringstream iss(line);
63 G4double prob;
64 iss >> prob;
65 sum_probs += prob;
66 probabilities.push_back(prob);
67 std::vector<std::string> types;
68 std::string type;
69 while (iss >> type) {
70 types.push_back(type);
71 }
72 particle_types.push_back(std::move(types));
73 }
74 }
75 else std::cout << "ERROR no fread_file " << filename << std::endl;
76 return sum_probs;
77}
78
79//this function will tell you the FS line number from the datafile based on your random value
80G4int NNbarToAnnihilationChannel::findStringNumber(G4double rdm, std::vector<G4double> yields) {
81 G4int stringNumber = -1;
82 G4double smallestsum = 0.0;
83 G4double biggestsum = yields[0];
84 //std::cout << "initial input " << rdm << std::endl;
85 for (G4int i = 0; i < static_cast<G4int>(yields.size()-1); i++) {
86 if (rdm >= smallestsum && rdm <= biggestsum) {
87 //std::cout << smallestsum << " and " << biggestsum << std::endl;
88 stringNumber = i+1;
89 }
90 smallestsum += yields[i];
91 biggestsum += yields[i+1];
92 }
93 if(stringNumber==-1) stringNumber = static_cast<G4int>(yields.size());
94 if(stringNumber==-1){
95 INCL_ERROR("ERROR in findStringNumber (stringNumber=-1)");
96 std::cout << "ERROR in findStringNumber" << std::endl;
97 }
98 return stringNumber;
99}
100
101
102void NNbarToAnnihilationChannel::fillFinalState(FinalState *fs) {
103
104 Particle *nucleon;
105 Particle *antinucleon;
106
107 if(particle1->isNucleon()){
108 nucleon = particle1;
109 antinucleon = particle2;
110 }
111 else{
112 nucleon = particle2;
113 antinucleon = particle1;
114 }
115
116 const G4double plab = 0.001*KinematicsUtils::momentumInLab(particle1, particle2); //GeV
117 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(nucleon, antinucleon);
118 const G4bool IsOnlyPion = (sqrtS < nucleon->getINCLMass() + antinucleon->getINCLMass());
119 G4double rdm = Random::shoot();
120
121 const std::vector<G4double> BFMM6 = {66.098, 0.153, -4.576, -38.319, 6.625}; //ppbar annihilation xs
122 const std::vector<G4double> BFMM1 = {119.066, 6.251, -0.006, -60.046, 11.958}; //ppbar total xs
123 const std::vector<G4double> BFMM471 = {108.104, 15.708, 0.832, -54.632, -6.958}; //npbar total xs
124 const std::vector<G4double> coef_nbarp_total = {1.69447, 5.26254E+08, -5.36346, -0.39766, 0.0243057}; //nbarp total xs (plab <0.5)
125
126 //PPbar annihilation xs
127 const std::vector<G4double> PPbar_pip_pim = {0.637, -0.340, -0.003, -0.439, 0.144};
128 const std::vector<G4double> PPbar_pip_pim_pi0 = {-2.065, 4.893, -1.130, 1.231, -0.212};
129 const std::vector<G4double> PPbar_pip_pim_omega = {3.020, 0.425, -0.029, -3.420, 0.867};
130 const std::vector<G4double> PPbar_pip_pim_Kp_Km = {-1.295, 1.897, -0.001, -0.365, 0.044};
131 const std::vector<G4double> PPbar_pip_pim_pi0_Kp_Km = {-12.220, 12.509, -0.351, 4.682, -0.777};
132 const std::vector<G4double> PPbar_2pip_2pim = {3.547, 0.095, 0.957, -3.444, 0.685};
133 const std::vector<G4double> PPbar_2pip_2pim_pi0 = {13.044, 1.449, 0.695, -12.313, 1.627};
134 const std::vector<G4double> PPbar_2pip_2pim_3pi0 = {6.398, 0.199, -1.103, -1.271, -0.380};
135 const std::vector<G4double> PPbar_3pip_3pim = {1.490, 0.240, 0.002, -1.012, 0.134};
136 const std::vector<G4double> PPbar_3pip_3pim_pi0 = {0.286, 1.634, -1.369, 3.099, -1.294};
137 const std::vector<G4double> PPbar_3pip_3pim_2pi0 = {-11.370, 12.503, -0.680, 10.059, -2.501};
138 const std::vector<G4double> PPbar_3pip_3pim_3pi0 = {-14.732, 12.338, -0.724, 11.342, -2.224};
139 const std::vector<G4double> PPbar_4pip_4pim = {-1.574, 1.607, -0.864, 1.253, -0.276};
140 const std::vector<G4double> PPbar_4pip_4pim_pi0 = {-1.096, 0.977, -0.995, 1.007, -0.171};
141
142 //NPbar annihilation xs
143 const std::vector<G4double> NPbar_pip_2pim = {-12.116, 14.485, -0.094, -1.632, 0.882, 5.000};
144 const std::vector<G4double> NPbar_pip_2pim_2pi0 = {8.276, 5.057, 0.483, -15.864, 2.552, 7.000};
145 const std::vector<G4double> NPbar_pip_2pim_3pi0 = {-1.500, 9.574, 0.528, -11.633, -0.615, 7.000};
146 const std::vector<G4double> NPbar_pip_2pim_pi0 = {7.999, 4.135, 0.608, -14.136, 1.590, 7.000};
147 const std::vector<G4double> NPbar_pip_pim_pi0_Km_K0 = {0.083, 0.091, -1.709, 0.284, -0.107};
148 const std::vector<G4double> NPbar_pip_pim_Km_K0 = {0.003, 0.297, -0.001, -0.143, 0.052};
149 const std::vector<G4double> NPbar_2pip_3pim_pi0 = {-14.701, 22.258, -0.001, -3.094, -0.190};
150 const std::vector<G4double> NPbar_2pip_3pim = {-0.616, 4.575, -0.002, -1.921, -0.153};
151
152
153 // File names
154 #ifdef INCLXX_IN_GEANT4_MODE
155 if(!G4FindDataDir("G4INCLDATA")) {
156 G4ExceptionDescription ed;
157 ed << " Data missing: set environment variable G4INCLDATA\n"
158 << " to point to the directory containing data files needed\n"
159 << " by the INCL++ model" << G4endl;
160 G4Exception("G4INCLDataFile::readData()","inflightppbarFS.dat, ...",
161 FatalException, ed);
162 }
163 G4String dataPath0{G4FindDataDir("G4INCLDATA")};
164 G4String dataPathppbar(dataPath0 + "/inflightppbarFS.dat");
165 G4String dataPathnpbar(dataPath0 + "/inflightnpbarFS.dat");
166 G4String dataPathppbark(dataPath0 + "/inflightppbarFSkaonic.dat");
167 G4String dataPathnpbark(dataPath0 + "/inflightnpbarFSkaonic.dat");
168 G4String dataPathnbarp(dataPath0 + "/inflightpnbarFS.dat"); //nbar case
169 G4String dataPathnbarpk(dataPath0 + "/inflightpnbarFSkaonic.dat"); // nbar case
170 #else
171 //Config *theConfig = new G4INCL::Config;
172 //theConfig->setINCLXXDataFilePath(G4INCL::theINCLXXDataFilePath);
173 Config const *theConfig=theNucleus->getStore()->getConfig();
174 std::string path;
175 if(theConfig)
176 path = theConfig->getINCLXXDataFilePath();
177 std::string dataPathppbar(path + "/inflightppbarFS.dat");
178 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
179 std::string dataPathnpbar(path + "/inflightnpbarFS.dat");
180 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
181 std::string dataPathppbark(path + "/inflightppbarFSkaonic.dat");
182 INCL_DEBUG("Reading https://doi.org/10.1016/j.physrep.2005.03.002 ppbar kaonic final states" << dataPathppbark << '\n');
183 std::string dataPathnpbark(path + "/inflightnpbarFSkaonic.dat");
184 INCL_DEBUG("Reading https://doi.org/10.1007/BF02818764 and https://link.springer.com/article/10.1007/BF02754930 npbar kaonic final states" << dataPathnpbark << '\n');
185 std::string dataPathnbarp(path + "/inflightpnbarFS.dat"); // nbar case
186 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N nnbar final states " << dataPathnbarp << '\n');
187 std::string dataPathnbarpk(path + "/inflightpnbarFSkaonic.dat");
188 INCL_DEBUG("Reading https://doi.org/10.1007/BF02818764 and https://link.springer.com/article/10.1007/BF02754930 nbarp kaonic final states" << dataPathnbarpk << '\n');
189 #endif
190 /*std::string path = {"/home/zdemid/INCL/inclcode/data"};
191 std::string dataPathppbar(path + "/inflightppbarFS.dat");
192 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
193 std::string dataPathnpbar(path + "/inflightnpbarFS.dat");
194 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
195 std::string dataPathppbark(path + "/inflightppbarFSkaonic.dat");
196 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar kaonic final states" << dataPathppbark << '\n');
197 std::string dataPathnpbark(path + "/inflightnpbarFSkaonic.dat");
198 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar kaonic final states" << dataPathnpbark << '\n');
199 every time we remove lines for which we have data from BFMM
200 */
201
202 std::vector<G4double> probabilities; //will store each FS yield
203 std::vector<std::vector<std::string>> particle_types; //will store particle names
204 G4double sum; //will contain a sum of probabilities of all FS in the file
205 const G4double kaonicFSprob=0.05; //probability to kave kaonic FS
206
207
208 ParticleList list;
209 //list.push_back(nucleon);
210 //list.push_back(antinucleon);
211 // NNbar will not be in the list because they annihilate
212 const ThreeVector &rcol = nucleon->getPosition();
213 const ThreeVector zero;
214
215 //setting types of new particles and pushing them back to the list
216 if(nucleon->getType()==Neutron && antinucleon->getType()==antiProton){
217 if(IsOnlyPion){
218 const std::vector<double> channels_Ratio = {0.053, 0.258, 0.531, 0.067, 0.062, 0.007, 0.020, 0.002};//Taken from inflightnpbarFS.dat and renormalised
219 if(rdm < channels_Ratio.front()){
220 Particle *p1 = new Particle(PiMinus, zero, rcol);
221 Particle *p2 = new Particle(PiZero, zero, rcol);
222
223 list.push_back(p1);
224 list.push_back(p2);
225 }
226 else if(rdm < channels_Ratio.front() + channels_Ratio[1]){
227 Particle *p1 = new Particle(PiMinus, zero, rcol);
228 Particle *p2 = new Particle(PiZero, zero, rcol);
229 Particle *p3 = new Particle(PiZero, zero, rcol);
230
231 list.push_back(p1);
232 list.push_back(p2);
233 list.push_back(p3);
234 }
235 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2]){
236 Particle *p1 = new Particle(PiMinus, zero, rcol);
237 Particle *p2 = new Particle(PiZero, zero, rcol);
238 Particle *p3 = new Particle(PiZero, zero, rcol);
239 Particle *p4 = new Particle(PiZero, zero, rcol);
240
241 list.push_back(p1);
242 list.push_back(p2);
243 list.push_back(p3);
244 list.push_back(p4);
245 }
246 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3]){
247 Particle *p1 = new Particle(PiMinus, zero, rcol);
248 Particle *p2 = new Particle(PiZero, zero, rcol);
249 Particle *p3 = new Particle(PiZero, zero, rcol);
250 Particle *p4 = new Particle(PiZero, zero, rcol);
251 Particle *p5 = new Particle(PiZero, zero, rcol);
252
253 list.push_back(p1);
254 list.push_back(p2);
255 list.push_back(p3);
256 list.push_back(p4);
257 list.push_back(p5);
258 }
259 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4]){
260 Particle *p1 = new Particle(PiMinus, zero, rcol);
261 Particle *p2 = new Particle(PiZero, zero, rcol);
262 Particle *p3 = new Particle(PiZero, zero, rcol);
263 Particle *p4 = new Particle(PiZero, zero, rcol);
264 Particle *p5 = new Particle(PiZero, zero, rcol);
265 Particle *p6 = new Particle(PiZero, zero, rcol);
266
267 list.push_back(p1);
268 list.push_back(p2);
269 list.push_back(p3);
270 list.push_back(p4);
271 list.push_back(p5);
272 list.push_back(p6);
273 }
274 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5]){
275 Particle *p1 = new Particle(PiPlus, zero, rcol);
276 Particle *p2 = new Particle(PiMinus, zero, rcol);
277 Particle *p3 = new Particle(PiMinus, zero, rcol);
278 Particle *p4 = new Particle(PiZero, zero, rcol);
279 Particle *p5 = new Particle(PiZero, zero, rcol);
280 Particle *p6 = new Particle(PiZero, zero, rcol);
281 Particle *p7 = new Particle(PiZero, zero, rcol);
282
283 list.push_back(p1);
284 list.push_back(p2);
285 list.push_back(p3);
286 list.push_back(p4);
287 list.push_back(p5);
288 list.push_back(p6);
289 list.push_back(p7);
290 }
291 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6]){
292 Particle *p1 = new Particle(PiPlus, zero, rcol);
293 Particle *p2 = new Particle(PiPlus, zero, rcol);
294 Particle *p3 = new Particle(PiMinus, zero, rcol);
295 Particle *p4 = new Particle(PiMinus, zero, rcol);
296 Particle *p5 = new Particle(PiMinus, zero, rcol);
297 Particle *p6 = new Particle(PiZero, zero, rcol);
298 Particle *p7 = new Particle(PiZero, zero, rcol);
299
300 list.push_back(p1);
301 list.push_back(p2);
302 list.push_back(p3);
303 list.push_back(p4);
304 list.push_back(p5);
305 list.push_back(p6);
306 list.push_back(p7);
307 }
308 else if(rdm <= channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6] + channels_Ratio[7]){
309 Particle *p1 = new Particle(PiPlus, zero, rcol);
310 Particle *p2 = new Particle(PiPlus, zero, rcol);
311 Particle *p3 = new Particle(PiPlus, zero, rcol);
312 Particle *p4 = new Particle(PiMinus, zero, rcol);
313 Particle *p5 = new Particle(PiMinus, zero, rcol);
314 Particle *p6 = new Particle(PiMinus, zero, rcol);
315 Particle *p7 = new Particle(PiMinus, zero, rcol);
316
317 list.push_back(p1);
318 list.push_back(p2);
319 list.push_back(p3);
320 list.push_back(p4);
321 list.push_back(p5);
322 list.push_back(p6);
323 list.push_back(p7);
324 }
325 else{
326 INCL_ERROR("random draw outside channels for OnlyPion annihilation (low energy)");
327 }
328 }
329 else {
330 const G4double totalpnbar = KinematicsUtils::compute_xs(BFMM6, plab)*KinematicsUtils::compute_xs(BFMM471, plab)/KinematicsUtils::compute_xs(BFMM1, plab);
331 // xs is same for npbar, but the fs has different charge
332
333 if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab)) {
334 // First condition
335 Particle* p1 = new Particle(PiPlus, zero, rcol);
336 Particle* p2 = new Particle(PiMinus, zero, rcol);
337 Particle* p3 = new Particle(PiMinus, zero, rcol);
338
339 list.push_back(p1);
340 list.push_back(p2);
341 list.push_back(p3);
342 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
343 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab)) {
344 // Second condition
345 Particle* p1 = new Particle(PiPlus, zero, rcol);
346 Particle* p2 = new Particle(PiMinus, zero, rcol);
347 Particle* p3 = new Particle(PiZero, zero, rcol);
348 Particle* p4 = new Particle(PiZero, zero, rcol);
349 Particle* p5 = new Particle(PiMinus, zero, rcol);
350
351 list.push_back(p1);
352 list.push_back(p2);
353 list.push_back(p3);
354 list.push_back(p4);
355 list.push_back(p5);
356 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
357 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
358 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab)) {
359 // Third condition
360 Particle* p1 = new Particle(PiPlus, zero, rcol);
361 Particle* p2 = new Particle(PiMinus, zero, rcol);
362 Particle* p3 = new Particle(PiZero, zero, rcol);
363 Particle* p4 = new Particle(PiZero, zero, rcol);
364 Particle* p5 = new Particle(PiZero, zero, rcol);
365 Particle* p6 = new Particle(PiMinus, zero, rcol);
366
367 list.push_back(p1);
368 list.push_back(p2);
369 list.push_back(p3);
370 list.push_back(p4);
371 list.push_back(p5);
372 list.push_back(p6);
373 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
374 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
375 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
376 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab)) {
377 // Fourth condition
378 Particle* p1 = new Particle(PiPlus, zero, rcol);
379 Particle* p2 = new Particle(PiMinus, zero, rcol);
380 Particle* p3 = new Particle(PiZero, zero, rcol);
381 Particle* p4 = new Particle(PiMinus, zero, rcol);
382
383 list.push_back(p1);
384 list.push_back(p2);
385 list.push_back(p3);
386 list.push_back(p4);
387 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
388 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
389 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
390 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
391 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab)) {
392 // Fifth condition
393 Particle* p1 = new Particle(PiPlus, zero, rcol);
394 Particle* p2 = new Particle(PiMinus, zero, rcol);
395 Particle* p3 = new Particle(PiZero, zero, rcol);
396 Particle* p4 = new Particle(KMinus, zero, rcol);
397 Particle* p5 = new Particle(KZero, zero, rcol);
398
399 list.push_back(p1);
400 list.push_back(p2);
401 list.push_back(p3);
402 list.push_back(p4);
403 list.push_back(p5);
404 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
405 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
406 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
407 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
408 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
409 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab)) {
410 // Sixth condition
411 Particle* p1 = new Particle(PiPlus, zero, rcol);
412 Particle* p2 = new Particle(PiMinus, zero, rcol);
413 Particle* p3 = new Particle(KMinus, zero, rcol);
414 Particle* p4 = new Particle(KZero, zero, rcol);
415
416 list.push_back(p1);
417 list.push_back(p2);
418 list.push_back(p3);
419 list.push_back(p4);
420 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
421 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
422 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
423 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
424 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
425 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
426 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab)) {
427 // Seventh condition
428 Particle* p1 = new Particle(PiPlus, zero, rcol);
429 Particle* p2 = new Particle(PiPlus, zero, rcol);
430 Particle* p3 = new Particle(PiMinus, zero, rcol);
431 Particle* p4 = new Particle(PiMinus, zero, rcol);
432 Particle* p5 = new Particle(PiMinus, zero, rcol);
433 Particle* p6 = new Particle(PiZero, zero, rcol);
434
435 list.push_back(p1);
436 list.push_back(p2);
437 list.push_back(p3);
438 list.push_back(p4);
439 list.push_back(p5);
440 list.push_back(p6);
441 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
442 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
443 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
444 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
445 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
446 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
447 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab) +
448 KinematicsUtils::compute_xs(NPbar_2pip_3pim, plab)) {
449 // Eighth condition
450 Particle* p1 = new Particle(PiPlus, zero, rcol);
451 Particle* p2 = new Particle(PiPlus, zero, rcol);
452 Particle* p3 = new Particle(PiMinus, zero, rcol);
453 Particle* p4 = new Particle(PiMinus, zero, rcol);
454 Particle* p5 = new Particle(PiMinus, zero, rcol);
455
456 list.push_back(p1);
457 list.push_back(p2);
458 list.push_back(p3);
459 list.push_back(p4);
460 list.push_back(p5);
461 } else {
462 // Default condition
463 //std::cout << "default condition pnbar"<< std::endl;
464 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
465 INCL_DEBUG("pionic npbar final state chosen" << '\n');
466 sum = read_file(dataPathnpbar, probabilities, particle_types);
467 rdm = (rdm/(1.-kaonicFSprob))*sum; //normalize by the sum of probabilities in the file
468 //now get the line number in the file where the FS particles are stored:
469 G4int n = findStringNumber(rdm, probabilities)-1;
470 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
471 if(particle_types[n][j] == "pi0"){
472 Particle *p = new Particle(PiZero, zero, rcol);
473 list.push_back(p);
474 }
475 else if(particle_types[n][j] == "pi-"){
476 Particle *p = new Particle(PiMinus, zero, rcol);
477 list.push_back(p);
478 }
479 else if(particle_types[n][j] == "pi+"){
480 Particle *p = new Particle(PiPlus, zero, rcol);
481 list.push_back(p);
482 }
483 else if(particle_types[n][j] == "omega"){
484 Particle *p = new Particle(Omega, zero, rcol);
485 list.push_back(p);
486 }
487 else if(particle_types[n][j] == "eta"){
488 Particle *p = new Particle(Eta, zero, rcol);
489 list.push_back(p);
490 }
491 else{
492 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
493 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
494 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
495 }
496 }
497 }
498 } // end of pionic option
499 else{
500 INCL_DEBUG("kaonic npbar final state chosen" << '\n');
501 sum = read_file(dataPathnpbark, probabilities, particle_types);
502 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
503 //now get the line number in the file where the FS particles are stored:
504 G4int n = findStringNumber(rdm, probabilities)-1;
505 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
506 if(particle_types[n][j] == "pi0"){
507 Particle *p = new Particle(PiZero, zero, rcol);
508 list.push_back(p);
509 }
510 else if(particle_types[n][j] == "pi-"){
511 Particle *p = new Particle(PiMinus, zero, rcol);
512 list.push_back(p);
513 }
514 else if(particle_types[n][j] == "pi+"){
515 Particle *p = new Particle(PiPlus, zero, rcol);
516 list.push_back(p);
517 }
518 else if(particle_types[n][j] == "omega"){
519 Particle *p = new Particle(Omega, zero, rcol);
520 list.push_back(p);
521 }
522 else if(particle_types[n][j] == "eta"){
523 Particle *p = new Particle(Eta, zero, rcol);
524 list.push_back(p);
525 }
526 else if(particle_types[n][j] == "K-"){
527 Particle *p = new Particle(KMinus, zero, rcol);
528 list.push_back(p);
529 }
530 else if(particle_types[n][j] == "K+"){
531 Particle *p = new Particle(KPlus, zero, rcol);
532 list.push_back(p);
533 }
534 else if(particle_types[n][j] == "K0"){
535 Particle *p = new Particle(KZero, zero, rcol);
536 list.push_back(p);
537 }
538 else if(particle_types[n][j] == "K0b"){
539 Particle *p = new Particle(KZeroBar, zero, rcol);
540 list.push_back(p);
541 }
542 else{
543 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
544 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
545 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
546 }
547 }
548 }
549 } // end of kaonic option
550 } // end of default annihilation
551 } //end of else for IsOnlyPion
552 }
553 else if(nucleon->getType()==Proton && antinucleon->getType()==antiNeutron){
554 if(IsOnlyPion){
555 const std::vector<G4double> channels_Ratio = {0.053, 0.258, 0.531, 0.067, 0.062, 0.007, 0.020, 0.002};//Taken from inflightpnbarFS.dat and renormalised
556 if(rdm < channels_Ratio.front()){
557 Particle *p1 = new Particle(PiPlus, zero, rcol);
558 Particle *p2 = new Particle(PiZero, zero, rcol);
559
560 list.push_back(p1);
561 list.push_back(p2);
562 }
563 else if(rdm < channels_Ratio.front() + channels_Ratio[1]){
564 Particle *p1 = new Particle(PiPlus, zero, rcol);
565 Particle *p2 = new Particle(PiZero, zero, rcol);
566 Particle *p3 = new Particle(PiZero, zero, rcol);
567
568 list.push_back(p1);
569 list.push_back(p2);
570 list.push_back(p3);
571 }
572 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2]){
573 Particle *p1 = new Particle(PiPlus, zero, rcol);
574 Particle *p2 = new Particle(PiZero, zero, rcol);
575 Particle *p3 = new Particle(PiZero, zero, rcol);
576 Particle *p4 = new Particle(PiZero, zero, rcol);
577
578 list.push_back(p1);
579 list.push_back(p2);
580 list.push_back(p3);
581 list.push_back(p4);
582 }
583 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3]){
584 Particle *p1 = new Particle(PiPlus, zero, rcol);
585 Particle *p2 = new Particle(PiZero, zero, rcol);
586 Particle *p3 = new Particle(PiZero, zero, rcol);
587 Particle *p4 = new Particle(PiZero, zero, rcol);
588 Particle *p5 = new Particle(PiZero, zero, rcol);
589
590 list.push_back(p1);
591 list.push_back(p2);
592 list.push_back(p3);
593 list.push_back(p4);
594 list.push_back(p5);
595 }
596 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4]){
597 Particle *p1 = new Particle(PiPlus, zero, rcol);
598 Particle *p2 = new Particle(PiZero, zero, rcol);
599 Particle *p3 = new Particle(PiZero, zero, rcol);
600 Particle *p4 = new Particle(PiZero, zero, rcol);
601 Particle *p5 = new Particle(PiZero, zero, rcol);
602 Particle *p6 = new Particle(PiZero, zero, rcol);
603
604 list.push_back(p1);
605 list.push_back(p2);
606 list.push_back(p3);
607 list.push_back(p4);
608 list.push_back(p5);
609 list.push_back(p6);
610 }
611 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5]){
612 Particle *p1 = new Particle(PiPlus, zero, rcol);
613 Particle *p2 = new Particle(PiMinus, zero, rcol);
614 Particle *p3 = new Particle(PiPlus, zero, rcol);
615 Particle *p4 = new Particle(PiZero, zero, rcol);
616 Particle *p5 = new Particle(PiZero, zero, rcol);
617 Particle *p6 = new Particle(PiZero, zero, rcol);
618 Particle *p7 = new Particle(PiZero, zero, rcol);
619
620 list.push_back(p1);
621 list.push_back(p2);
622 list.push_back(p3);
623 list.push_back(p4);
624 list.push_back(p5);
625 list.push_back(p6);
626 list.push_back(p7);
627 }
628 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6]){
629 Particle *p1 = new Particle(PiPlus, zero, rcol);
630 Particle *p2 = new Particle(PiPlus, zero, rcol);
631 Particle *p3 = new Particle(PiMinus, zero, rcol);
632 Particle *p4 = new Particle(PiMinus, zero, rcol);
633 Particle *p5 = new Particle(PiPlus, zero, rcol);
634 Particle *p6 = new Particle(PiZero, zero, rcol);
635 Particle *p7 = new Particle(PiZero, zero, rcol);
636
637 list.push_back(p1);
638 list.push_back(p2);
639 list.push_back(p3);
640 list.push_back(p4);
641 list.push_back(p5);
642 list.push_back(p6);
643 list.push_back(p7);
644 }
645 else if(rdm <= channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6] + channels_Ratio[7]){
646 Particle *p1 = new Particle(PiPlus, zero, rcol);
647 Particle *p2 = new Particle(PiPlus, zero, rcol);
648 Particle *p3 = new Particle(PiPlus, zero, rcol);
649 Particle *p4 = new Particle(PiMinus, zero, rcol);
650 Particle *p5 = new Particle(PiMinus, zero, rcol);
651 Particle *p6 = new Particle(PiMinus, zero, rcol);
652 Particle *p7 = new Particle(PiPlus, zero, rcol);
653
654 list.push_back(p1);
655 list.push_back(p2);
656 list.push_back(p3);
657 list.push_back(p4);
658 list.push_back(p5);
659 list.push_back(p6);
660 list.push_back(p7);
661 }
662 }
663 else {
664 if (plab < 1.){ //nbar != pbar (< 1. GeV/c)
665 if(rdm < 1. - kaonicFSprob){
666 INCL_DEBUG("pionic pnbar final state" << '\n');
667 sum = read_file(dataPathnbarp, probabilities, particle_types);
668 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.95 normalize by the sum of probabilities in the file
669 //now get the line number in the file where the FS particles are stored:
670 G4int n = findStringNumber(rdm, probabilities)-1;
671 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
672 if(particle_types[n][j] == "pi0"){
673 Particle *p = new Particle(PiZero, zero, rcol);
674 list.push_back(p);
675 }
676 else if(particle_types[n][j] == "pi-"){
677 Particle *p = new Particle(PiMinus, zero, rcol);
678 list.push_back(p);
679 }
680 else if(particle_types[n][j] == "pi+"){
681 Particle *p = new Particle(PiPlus, zero, rcol);
682 list.push_back(p);
683 }
684 else if(particle_types[n][j] == "omega"){
685 Particle *p = new Particle(Omega, zero, rcol);
686 list.push_back(p);
687 }
688 else if(particle_types[n][j] == "eta"){
689 Particle *p = new Particle(Eta, zero, rcol);
690 list.push_back(p);
691 }
692 else{
693 INCL_ERROR("Some non-existing FS particle detected when reading nbar FS files");
694 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
695 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
696 }
697 }
698 }
699 }
700 else{
701 INCL_DEBUG("kaonic npbar final state chosen" << '\n');
702 sum = read_file(dataPathnbarpk, probabilities, particle_types);
703 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
704 //now get the line number in the file where the FS particles are stored:
705 G4int n = findStringNumber(rdm, probabilities)-1;
706 for(G4int j = 0; j < static_cast<int>(particle_types[n].size()); j++){
707 if(particle_types[n][j] == "pi0"){
708 Particle *p = new Particle(PiZero, zero, rcol);
709 list.push_back(p);
710 }
711 else if(particle_types[n][j] == "pi-"){
712 Particle *p = new Particle(PiMinus, zero, rcol);
713 list.push_back(p);
714 }
715 else if(particle_types[n][j] == "pi+"){
716 Particle *p = new Particle(PiPlus, zero, rcol);
717 list.push_back(p);
718 }
719 else if(particle_types[n][j] == "omega"){
720 Particle *p = new Particle(Omega, zero, rcol);
721 list.push_back(p);
722 }
723 else if(particle_types[n][j] == "eta"){
724 Particle *p = new Particle(Eta, zero, rcol);
725 list.push_back(p);
726 }
727 else if(particle_types[n][j] == "K-"){
728 Particle *p = new Particle(KMinus, zero, rcol);
729 list.push_back(p);
730 }
731 else if(particle_types[n][j] == "K+"){
732 Particle *p = new Particle(KPlus, zero, rcol);
733 list.push_back(p);
734 }
735 else if(particle_types[n][j] == "K0"){
736 Particle *p = new Particle(KZero, zero, rcol);
737 list.push_back(p);
738 }
739 else if(particle_types[n][j] == "K0b"){
740 Particle *p = new Particle(KZeroBar, zero, rcol);
741 list.push_back(p);
742 }
743 else{
744 INCL_ERROR("Some non-existing FS particle detected when reading nbar FS files");
745 for(G4int jj = 0; jj < static_cast<int>(particle_types[n].size()); jj++){
746 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
747 }
748 }
749 }
750 }
751 }
752 else{ //nbar = pbar (>1000 MeV/c)
753 const G4double totalpnbar = KinematicsUtils::compute_xs(BFMM6, plab)*KinematicsUtils::compute_xs(BFMM471, plab)/KinematicsUtils::compute_xs(BFMM1, plab);
754 // xs is same for npbar, but the fs has different charge
755
756 if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab)) {
757 // First condition
758 Particle* p1 = new Particle(PiPlus, zero, rcol);
759 Particle* p2 = new Particle(PiMinus, zero, rcol);
760 Particle* p3 = new Particle(PiPlus, zero, rcol);
761
762 list.push_back(p1);
763 list.push_back(p2);
764 list.push_back(p3);
765 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
766 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab)) {
767 // Second condition
768 Particle* p1 = new Particle(PiPlus, zero, rcol);
769 Particle* p2 = new Particle(PiMinus, zero, rcol);
770 Particle* p3 = new Particle(PiZero, zero, rcol);
771 Particle* p4 = new Particle(PiZero, zero, rcol);
772 Particle* p5 = new Particle(PiPlus, zero, rcol);
773
774 list.push_back(p1);
775 list.push_back(p2);
776 list.push_back(p3);
777 list.push_back(p4);
778 list.push_back(p5);
779 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
780 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
781 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab)) {
782 // Third condition
783 Particle* p1 = new Particle(PiPlus, zero, rcol);
784 Particle* p2 = new Particle(PiMinus, zero, rcol);
785 Particle* p3 = new Particle(PiZero, zero, rcol);
786 Particle* p4 = new Particle(PiZero, zero, rcol);
787 Particle* p5 = new Particle(PiZero, zero, rcol);
788 Particle* p6 = new Particle(PiPlus, zero, rcol);
789
790 list.push_back(p1);
791 list.push_back(p2);
792 list.push_back(p3);
793 list.push_back(p4);
794 list.push_back(p5);
795 list.push_back(p6);
796 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
797 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
798 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
799 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab)) {
800 // Fourth condition
801 Particle* p1 = new Particle(PiPlus, zero, rcol);
802 Particle* p2 = new Particle(PiMinus, zero, rcol);
803 Particle* p3 = new Particle(PiZero, zero, rcol);
804 Particle* p4 = new Particle(PiPlus, zero, rcol);
805
806 list.push_back(p1);
807 list.push_back(p2);
808 list.push_back(p3);
809 list.push_back(p4);
810 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
811 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
812 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
813 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
814 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab)) {
815 // Fifth condition
816 Particle* p1 = new Particle(PiPlus, zero, rcol);
817 Particle* p2 = new Particle(PiMinus, zero, rcol);
818 Particle* p3 = new Particle(PiZero, zero, rcol);
819 Particle* p4 = new Particle(KPlus, zero, rcol);
820 Particle* p5 = new Particle(KZeroBar, zero, rcol);
821
822 list.push_back(p1);
823 list.push_back(p2);
824 list.push_back(p3);
825 list.push_back(p4);
826 list.push_back(p5);
827 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
828 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
829 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
830 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
831 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
832 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab)) {
833 // Sixth condition
834 Particle* p1 = new Particle(PiPlus, zero, rcol);
835 Particle* p2 = new Particle(PiMinus, zero, rcol);
836 Particle* p3 = new Particle(KPlus, zero, rcol);
837 Particle* p4 = new Particle(KZeroBar, zero, rcol);
838
839 list.push_back(p1);
840 list.push_back(p2);
841 list.push_back(p3);
842 list.push_back(p4);
843 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
844 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
845 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
846 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
847 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
848 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
849 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab)) {
850 // Seventh condition
851 Particle* p1 = new Particle(PiPlus, zero, rcol);
852 Particle* p2 = new Particle(PiPlus, zero, rcol);
853 Particle* p3 = new Particle(PiMinus, zero, rcol);
854 Particle* p4 = new Particle(PiMinus, zero, rcol);
855 Particle* p5 = new Particle(PiPlus, zero, rcol);
856 Particle* p6 = new Particle(PiZero, zero, rcol);
857
858 list.push_back(p1);
859 list.push_back(p2);
860 list.push_back(p3);
861 list.push_back(p4);
862 list.push_back(p5);
863 list.push_back(p6);
864 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
865 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
866 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
867 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
868 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
869 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
870 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab) +
871 KinematicsUtils::compute_xs(NPbar_2pip_3pim, plab)) {
872 // Eighth condition
873 Particle* p1 = new Particle(PiPlus, zero, rcol);
874 Particle* p2 = new Particle(PiPlus, zero, rcol);
875 Particle* p3 = new Particle(PiMinus, zero, rcol);
876 Particle* p4 = new Particle(PiMinus, zero, rcol);
877 Particle* p5 = new Particle(PiPlus, zero, rcol);
878
879 list.push_back(p1);
880 list.push_back(p2);
881 list.push_back(p3);
882 list.push_back(p4);
883 list.push_back(p5);
884 } else {
885 // Default condition
886 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
887 INCL_DEBUG("pionic pnbar final state chosen" << '\n');
888 sum = read_file(dataPathnbarp, probabilities, particle_types);
889 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.95 normalize by the sum of probabilities in the file
890 //now get the line number in the file where the FS particles are stored:
891 G4int n = findStringNumber(rdm, probabilities)-1;
892 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
893 if(particle_types[n][j] == "pi0"){
894 Particle *p = new Particle(PiZero, zero, rcol);
895 list.push_back(p);
896 }
897 else if(particle_types[n][j] == "pi-"){
898 Particle *p = new Particle(PiMinus, zero, rcol);
899 list.push_back(p);
900 }
901 else if(particle_types[n][j] == "pi+"){
902 Particle *p = new Particle(PiPlus, zero, rcol);
903 list.push_back(p);
904 }
905 else if(particle_types[n][j] == "omega"){
906 Particle *p = new Particle(Omega, zero, rcol);
907 list.push_back(p);
908 }
909 else if(particle_types[n][j] == "eta"){
910 Particle *p = new Particle(Eta, zero, rcol);
911 list.push_back(p);
912 }
913 else{
914 INCL_ERROR("Some non-existing FS particle detected when reading nbar FS files");
915 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
916 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
917 }
918 }
919 }
920 } // end of pionic option
921 else{
922 INCL_DEBUG("kaonic pnbar final state chosen" << '\n');
923 sum = read_file(dataPathnbarpk, probabilities, particle_types);
924 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
925 //now get the line number in the file where the FS particles are stored:
926 G4int n = findStringNumber(rdm, probabilities)-1;
927 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
928 if(particle_types[n][j] == "pi0"){
929 Particle *p = new Particle(PiZero, zero, rcol);
930 list.push_back(p);
931 }
932 else if(particle_types[n][j] == "pi-"){
933 Particle *p = new Particle(PiMinus, zero, rcol);
934 list.push_back(p);
935 }
936 else if(particle_types[n][j] == "pi+"){
937 Particle *p = new Particle(PiPlus, zero, rcol);
938 list.push_back(p);
939 }
940 else if(particle_types[n][j] == "omega"){
941 Particle *p = new Particle(Omega, zero, rcol);
942 list.push_back(p);
943 }
944 else if(particle_types[n][j] == "eta"){
945 Particle *p = new Particle(Eta, zero, rcol);
946 list.push_back(p);
947 }
948 else if(particle_types[n][j] == "K-"){
949 Particle *p = new Particle(KMinus, zero, rcol);
950 list.push_back(p);
951 }
952 else if(particle_types[n][j] == "K+"){
953 Particle *p = new Particle(KPlus, zero, rcol);
954 list.push_back(p);
955 }
956 else if(particle_types[n][j] == "K0"){
957 Particle *p = new Particle(KZero, zero, rcol);
958 list.push_back(p);
959 }
960 else if(particle_types[n][j] == "K0b"){
961 Particle *p = new Particle(KZeroBar, zero, rcol);
962 list.push_back(p);
963 }
964 else{
965 INCL_ERROR("Some non-existing FS particle detected when reading pnbar FS files");
966 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
967 std::cout << "gotcha! " << particle_types[n][jj] << std::endl; }
968 }
969 }
970 } // end of kaonic option
971 } // end of default annihilation
972 }//end of nbarp & pbarn ( with Plab < 1 GeV/c)
973 } //end of else for IsOnlyPion
974 }
975 else{ //ppbar or nnbar
976 if(IsOnlyPion){
977 const std::vector<G4double> channels_Ratio = {0.0005, 0.0278, 0.0052, 0.3058, 0.0017, 0.3346, 0.0017, 0.0688, 0.2534, 0.0005};//Taken from inflightppbarFS.dat and renormalised
978 if(rdm < channels_Ratio.front()){
979 Particle *p1 = new Particle(PiZero, zero, rcol);
980 Particle *p2 = new Particle(PiZero, zero, rcol);
981
982 list.push_back(p1);
983 list.push_back(p2);
984 }
985 else if(rdm < channels_Ratio.front() + channels_Ratio[1]){
986 Particle *p1 = new Particle(PiZero, zero, rcol);
987 Particle *p2 = new Particle(PiZero, zero, rcol);
988 Particle *p3 = new Particle(PiZero, zero, rcol);
989
990 list.push_back(p1);
991 list.push_back(p2);
992 list.push_back(p3);
993 }
994 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2]){
995 Particle *p1 = new Particle(PiZero, zero, rcol);
996 Particle *p2 = new Particle(PiZero, zero, rcol);
997 Particle *p3 = new Particle(PiZero, zero, rcol);
998 Particle *p4 = new Particle(PiZero, zero, rcol);
999
1000 list.push_back(p1);
1001 list.push_back(p2);
1002 list.push_back(p3);
1003 list.push_back(p4);
1004 }
1005 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] ){
1006 Particle *p1 = new Particle(PiPlus, zero, rcol);
1007 Particle *p2 = new Particle(PiMinus, zero, rcol);
1008 Particle *p3 = new Particle(PiZero, zero, rcol);
1009 Particle *p4 = new Particle(PiZero, zero, rcol);
1010
1011 list.push_back(p1);
1012 list.push_back(p2);
1013 list.push_back(p3);
1014 list.push_back(p4);
1015 }
1016 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4]){
1017 Particle *p1 = new Particle(PiZero, zero, rcol);
1018 Particle *p2 = new Particle(PiZero, zero, rcol);
1019 Particle *p3 = new Particle(PiZero, zero, rcol);
1020 Particle *p4 = new Particle(PiZero, zero, rcol);
1021 Particle *p5 = new Particle(PiZero, zero, rcol);
1022
1023 list.push_back(p1);
1024 list.push_back(p2);
1025 list.push_back(p3);
1026 list.push_back(p4);
1027 list.push_back(p5);
1028 }
1029 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5]){
1030 Particle *p1 = new Particle(PiPlus, zero, rcol);
1031 Particle *p2 = new Particle(PiMinus, zero, rcol);
1032 Particle *p3 = new Particle(PiZero, zero, rcol);
1033 Particle *p4 = new Particle(PiZero, zero, rcol);
1034 Particle *p5 = new Particle(PiZero, zero, rcol);
1035
1036 list.push_back(p1);
1037 list.push_back(p2);
1038 list.push_back(p3);
1039 list.push_back(p4);
1040 list.push_back(p5);
1041 }
1042 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6]){
1043 Particle *p1 = new Particle(PiZero, zero, rcol);
1044 Particle *p2 = new Particle(PiZero, zero, rcol);
1045 Particle *p3 = new Particle(PiZero, zero, rcol);
1046 Particle *p4 = new Particle(PiZero, zero, rcol);
1047 Particle *p5 = new Particle(PiZero, zero, rcol);
1048 Particle *p6 = new Particle(PiZero, zero, rcol);
1049
1050 list.push_back(p1);
1051 list.push_back(p2);
1052 list.push_back(p3);
1053 list.push_back(p4);
1054 list.push_back(p5);
1055 list.push_back(p6);
1056 }
1057 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6] + channels_Ratio[7]){
1058 Particle *p1 = new Particle(PiPlus, zero, rcol);
1059 Particle *p2 = new Particle(PiMinus, zero, rcol);
1060 Particle *p3 = new Particle(PiZero, zero, rcol);
1061 Particle *p4 = new Particle(PiZero, zero, rcol);
1062 Particle *p5 = new Particle(PiZero, zero, rcol);
1063 Particle *p6 = new Particle(PiZero, zero, rcol);
1064
1065 list.push_back(p1);
1066 list.push_back(p2);
1067 list.push_back(p3);
1068 list.push_back(p4);
1069 list.push_back(p5);
1070 list.push_back(p6);
1071 }
1072 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6] + channels_Ratio[7] + channels_Ratio[8]){
1073 Particle *p1 = new Particle(PiPlus, zero, rcol);
1074 Particle *p2 = new Particle(PiPlus, zero, rcol);
1075 Particle *p3 = new Particle(PiMinus, zero, rcol);
1076 Particle *p4 = new Particle(PiMinus, zero, rcol);
1077 Particle *p5 = new Particle(PiZero, zero, rcol);
1078 Particle *p6 = new Particle(PiZero, zero, rcol);
1079
1080 list.push_back(p1);
1081 list.push_back(p2);
1082 list.push_back(p3);
1083 list.push_back(p4);
1084 list.push_back(p5);
1085 list.push_back(p6);
1086 }
1087 else if(rdm < channels_Ratio.front() + channels_Ratio[1] + channels_Ratio[2] + channels_Ratio[3] + channels_Ratio[4] + channels_Ratio[5] + channels_Ratio[6] + channels_Ratio[7] + channels_Ratio[8] + channels_Ratio[9]){
1088 Particle *p1 = new Particle(PiPlus, zero, rcol);
1089 Particle *p2 = new Particle(PiMinus, zero, rcol);
1090 Particle *p3 = new Particle(PiZero, zero, rcol);
1091 Particle *p4 = new Particle(PiZero, zero, rcol);
1092 Particle *p5 = new Particle(PiZero, zero, rcol);
1093 Particle *p6 = new Particle(PiZero, zero, rcol);
1094 Particle *p7 = new Particle(PiZero, zero, rcol);
1095
1096 list.push_back(p1);
1097 list.push_back(p2);
1098 list.push_back(p3);
1099 list.push_back(p4);
1100 list.push_back(p5);
1101 list.push_back(p6);
1102 list.push_back(p7);
1103 }
1104 else{
1105 INCL_ERROR("random draw outside channels for OnlyPion annihilation (low energy)");
1106 }
1107 }
1108 else{
1109 const G4double totalppbar = KinematicsUtils::compute_xs(BFMM6, plab);
1110 // same for nnbar
1111
1112 if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab)) {
1113 // First condition
1114 Particle* p1 = new Particle(PiPlus, zero, rcol);
1115 Particle* p2 = new Particle(PiMinus, zero, rcol);
1116
1117 list.push_back(p1);
1118 list.push_back(p2);
1119
1120
1121 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1122 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab)) {
1123 // Second condition
1124 Particle* p1 = new Particle(PiPlus, zero, rcol);
1125 Particle* p2 = new Particle(PiMinus, zero, rcol);
1126 Particle* p3 = new Particle(PiZero, zero, rcol);
1127
1128 list.push_back(p1);
1129 list.push_back(p2);
1130 list.push_back(p3);
1131 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1132 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1133 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab)) {
1134 // Third condition
1135 Particle* p1 = new Particle(PiPlus, zero, rcol);
1136 Particle* p2 = new Particle(PiMinus, zero, rcol);
1137 Particle* p3 = new Particle(Omega, zero, rcol);
1138
1139 list.push_back(p1);
1140 list.push_back(p2);
1141 list.push_back(p3);
1142 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1143 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1144 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1145 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab)) {
1146 // Fourth condition
1147 Particle* p1 = new Particle(PiPlus, zero, rcol);
1148 Particle* p2 = new Particle(PiMinus, zero, rcol);
1149 Particle* p3 = new Particle(KPlus, zero, rcol);
1150 Particle* p4 = new Particle(KMinus, zero, rcol);
1151
1152 list.push_back(p1);
1153 list.push_back(p2);
1154 list.push_back(p3);
1155 list.push_back(p4);
1156 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1157 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1158 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1159 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1160 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab)) {
1161 // Fifth condition
1162 Particle* p1 = new Particle(PiPlus, zero, rcol);
1163 Particle* p2 = new Particle(PiMinus, zero, rcol);
1164 Particle* p3 = new Particle(PiZero, zero, rcol);
1165 Particle* p4 = new Particle(KPlus, zero, rcol);
1166 Particle* p5 = new Particle(KMinus, zero, rcol);
1167
1168 list.push_back(p1);
1169 list.push_back(p2);
1170 list.push_back(p3);
1171 list.push_back(p4);
1172 list.push_back(p5);
1173 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1174 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1175 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1176 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1177 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1178 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab)) {
1179 // Sixth condition
1180 Particle* p1 = new Particle(PiPlus, zero, rcol);
1181 Particle* p2 = new Particle(PiMinus, zero, rcol);
1182 Particle* p3 = new Particle(PiPlus, zero, rcol);
1183 Particle* p4 = new Particle(PiMinus, zero, rcol);
1184
1185 list.push_back(p1);
1186 list.push_back(p2);
1187 list.push_back(p3);
1188 list.push_back(p4);
1189 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1190 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1191 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1192 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1193 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1194 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1195 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab)) {
1196 // Seventh condition
1197 Particle* p1 = new Particle(PiPlus, zero, rcol);
1198 Particle* p2 = new Particle(PiMinus, zero, rcol);
1199 Particle* p3 = new Particle(PiPlus, zero, rcol);
1200 Particle* p4 = new Particle(PiMinus, zero, rcol);
1201 Particle* p5 = new Particle(PiZero, zero, rcol);
1202
1203 list.push_back(p1);
1204 list.push_back(p2);
1205 list.push_back(p3);
1206 list.push_back(p4);
1207 list.push_back(p5);
1208 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1209 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1210 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1211 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1212 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1213 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1214 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
1215 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab)) {
1216 // Eighth condition
1217 Particle* p1 = new Particle(PiPlus, zero, rcol);
1218 Particle* p2 = new Particle(PiMinus, zero, rcol);
1219 Particle* p3 = new Particle(PiPlus, zero, rcol);
1220 Particle* p4 = new Particle(PiMinus, zero, rcol);
1221 Particle* p5 = new Particle(PiZero, zero, rcol);
1222 Particle* p6 = new Particle(PiZero, zero, rcol);
1223 Particle* p7 = new Particle(PiZero, zero, rcol);
1224
1225 list.push_back(p1);
1226 list.push_back(p2);
1227 list.push_back(p3);
1228 list.push_back(p4);
1229 list.push_back(p5);
1230 list.push_back(p6);
1231 list.push_back(p7);
1232 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1233 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1234 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1235 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1236 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1237 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1238 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
1239 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
1240 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab)) {
1241 // Ninth condition
1242 Particle* p1 = new Particle(PiPlus, zero, rcol);
1243 Particle* p2 = new Particle(PiMinus, zero, rcol);
1244 Particle* p3 = new Particle(PiPlus, zero, rcol);
1245 Particle* p4 = new Particle(PiMinus, zero, rcol);
1246 Particle* p5 = new Particle(PiPlus, zero, rcol);
1247 Particle* p6 = new Particle(PiMinus, zero, rcol);
1248
1249 list.push_back(p1);
1250 list.push_back(p2);
1251 list.push_back(p3);
1252 list.push_back(p4);
1253 list.push_back(p5);
1254 list.push_back(p6);
1255 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1256 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1257 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1258 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1259 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1260 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1261 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
1262 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
1263 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
1264 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab)) {
1265 // Tenth condition
1266 Particle* p1 = new Particle(PiPlus, zero, rcol);
1267 Particle* p2 = new Particle(PiMinus, zero, rcol);
1268 Particle* p3 = new Particle(PiPlus, zero, rcol);
1269 Particle* p4 = new Particle(PiMinus, zero, rcol);
1270 Particle* p5 = new Particle(PiPlus, zero, rcol);
1271 Particle* p6 = new Particle(PiMinus, zero, rcol);
1272 Particle* p7 = new Particle(PiZero, zero, rcol);
1273
1274 list.push_back(p1);
1275 list.push_back(p2);
1276 list.push_back(p3);
1277 list.push_back(p4);
1278 list.push_back(p5);
1279 list.push_back(p6);
1280 list.push_back(p7);
1281 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1282 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1283 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1284 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1285 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1286 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1287 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
1288 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
1289 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
1290 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
1291 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab)) {
1292 // Eleventh condition
1293 Particle* p1 = new Particle(PiPlus, zero, rcol);
1294 Particle* p2 = new Particle(PiMinus, zero, rcol);
1295 Particle* p3 = new Particle(PiPlus, zero, rcol);
1296 Particle* p4 = new Particle(PiMinus, zero, rcol);
1297 Particle* p5 = new Particle(PiPlus, zero, rcol);
1298 Particle* p6 = new Particle(PiMinus, zero, rcol);
1299 Particle* p7 = new Particle(PiZero, zero, rcol);
1300 Particle* p8 = new Particle(PiZero, zero, rcol);
1301
1302 list.push_back(p1);
1303 list.push_back(p2);
1304 list.push_back(p3);
1305 list.push_back(p4);
1306 list.push_back(p5);
1307 list.push_back(p6);
1308 list.push_back(p7);
1309 list.push_back(p8);
1310 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1311 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1312 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1313 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1314 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1315 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1316 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
1317 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
1318 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
1319 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
1320 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
1321 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab)) {
1322 // Twelfth condition
1323 Particle* p1 = new Particle(PiPlus, zero, rcol);
1324 Particle* p2 = new Particle(PiMinus, zero, rcol);
1325 Particle* p3 = new Particle(PiPlus, zero, rcol);
1326 Particle* p4 = new Particle(PiMinus, zero, rcol);
1327 Particle* p5 = new Particle(PiPlus, zero, rcol);
1328 Particle* p6 = new Particle(PiMinus, zero, rcol);
1329 Particle* p7 = new Particle(PiZero, zero, rcol);
1330 Particle* p8 = new Particle(PiZero, zero, rcol);
1331 Particle* p9 = new Particle(PiZero, zero, rcol);
1332
1333 list.push_back(p1);
1334 list.push_back(p2);
1335 list.push_back(p3);
1336 list.push_back(p4);
1337 list.push_back(p5);
1338 list.push_back(p6);
1339 list.push_back(p7);
1340 list.push_back(p8);
1341 list.push_back(p9);
1342 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1343 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1344 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1345 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1346 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1347 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1348 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
1349 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
1350 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
1351 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
1352 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
1353 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab) +
1354 KinematicsUtils::compute_xs(PPbar_4pip_4pim, plab)) {
1355 // Thirteenth condition
1356 Particle* p1 = new Particle(PiPlus, zero, rcol);
1357 Particle* p2 = new Particle(PiMinus, zero, rcol);
1358 Particle* p3 = new Particle(PiPlus, zero, rcol);
1359 Particle* p4 = new Particle(PiMinus, zero, rcol);
1360 Particle* p5 = new Particle(PiPlus, zero, rcol);
1361 Particle* p6 = new Particle(PiMinus, zero, rcol);
1362 Particle* p7 = new Particle(PiPlus, zero, rcol);
1363 Particle* p8 = new Particle(PiMinus, zero, rcol);
1364
1365 list.push_back(p1);
1366 list.push_back(p2);
1367 list.push_back(p3);
1368 list.push_back(p4);
1369 list.push_back(p5);
1370 list.push_back(p6);
1371 list.push_back(p7);
1372 list.push_back(p8);
1373 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
1374 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
1375 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
1376 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
1377 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
1378 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
1379 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
1380 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
1381 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
1382 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
1383 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
1384 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab) +
1385 KinematicsUtils::compute_xs(PPbar_4pip_4pim, plab) +
1386 KinematicsUtils::compute_xs(PPbar_4pip_4pim_pi0, plab)) {
1387 // Fourteenth condition
1388 Particle* p1 = new Particle(PiPlus, zero, rcol);
1389 Particle* p2 = new Particle(PiMinus, zero, rcol);
1390 Particle* p3 = new Particle(PiPlus, zero, rcol);
1391 Particle* p4 = new Particle(PiMinus, zero, rcol);
1392 Particle* p5 = new Particle(PiPlus, zero, rcol);
1393 Particle* p6 = new Particle(PiMinus, zero, rcol);
1394 Particle* p7 = new Particle(PiPlus, zero, rcol);
1395 Particle* p8 = new Particle(PiMinus, zero, rcol);
1396 Particle* p9 = new Particle(PiZero, zero, rcol);
1397
1398 list.push_back(p1);
1399 list.push_back(p2);
1400 list.push_back(p3);
1401 list.push_back(p4);
1402 list.push_back(p5);
1403 list.push_back(p6);
1404 list.push_back(p7);
1405 list.push_back(p8);
1406 list.push_back(p9);
1407 } else {
1408 // Default condition
1409 if(rdm < (1.-kaonicFSprob)){ // pionic FS was chosen
1410 INCL_DEBUG("pionic pp final state chosen" << '\n');
1411 sum = read_file(dataPathppbar, probabilities, particle_types);
1412 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.88 normalize by the sum of probabilities in the file
1413 //now get the line number in the file where the FS particles are stored:
1414 G4int n = findStringNumber(rdm, probabilities)-1;
1415 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
1416 if(particle_types[n][j] == "pi0"){
1417 Particle *p = new Particle(PiZero, zero, rcol);
1418 list.push_back(p);
1419 }
1420 else if(particle_types[n][j] == "pi-"){
1421 Particle *p = new Particle(PiMinus, zero, rcol);
1422 list.push_back(p);
1423 }
1424 else if(particle_types[n][j] == "pi+"){
1425 Particle *p = new Particle(PiPlus, zero, rcol);
1426 list.push_back(p);
1427 }
1428 else if(particle_types[n][j] == "omega"){
1429 Particle *p = new Particle(Omega, zero, rcol);
1430 list.push_back(p);
1431 }
1432 else if(particle_types[n][j] == "eta"){
1433 Particle *p = new Particle(Eta, zero, rcol);
1434 list.push_back(p);
1435 }
1436 else{
1437 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
1438 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
1439 std::cout << "gotcha! " << particle_types[n][jj] << std::endl; }
1440 }
1441 }
1442 } //end of pionic option
1443 else{
1444 INCL_DEBUG("kaonic pp final state chosen" << '\n');
1445 sum = read_file(dataPathppbark, probabilities, particle_types);
1446 rdm = ((1-rdm)/kaonicFSprob)*sum;//2670 normalize by the sum of probabilities in the file
1447 //now get the line number in the file where the FS particles are stored:
1448 G4int n = findStringNumber(rdm, probabilities)-1;
1449 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
1450 if(particle_types[n][j] == "pi0"){
1451 Particle *p = new Particle(PiZero, zero, rcol);
1452 list.push_back(p);
1453 }
1454 else if(particle_types[n][j] == "pi-"){
1455 Particle *p = new Particle(PiMinus, zero, rcol);
1456 list.push_back(p);
1457 }
1458 else if(particle_types[n][j] == "pi+"){
1459 Particle *p = new Particle(PiPlus, zero, rcol);
1460 list.push_back(p);
1461 }
1462 else if(particle_types[n][j] == "omega"){
1463 Particle *p = new Particle(Omega, zero, rcol);
1464 list.push_back(p);
1465 }
1466 else if(particle_types[n][j] == "eta"){
1467 Particle *p = new Particle(Eta, zero, rcol);
1468 list.push_back(p);
1469 }
1470 else if(particle_types[n][j] == "K-"){
1471 Particle *p = new Particle(KMinus, zero, rcol);
1472 list.push_back(p);
1473 }
1474 else if(particle_types[n][j] == "K+"){
1475 Particle *p = new Particle(KPlus, zero, rcol);
1476 list.push_back(p);
1477 }
1478 else if(particle_types[n][j] == "K0"){
1479 Particle *p = new Particle(KZero, zero, rcol);
1480 list.push_back(p);
1481 }
1482 else if(particle_types[n][j] == "K0b"){
1483 Particle *p = new Particle(KZeroBar, zero, rcol);
1484 list.push_back(p);
1485 }
1486 else{
1487 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
1488 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
1489 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
1490 }
1491 }
1492 }
1493 } // end of kaonic option
1494 } // end of default condition
1495 } //end of else for IsOnlyPion
1496 } // end of ppbar and nnbar case
1497
1498
1499
1500 nucleon->setType(list[0]->getType());
1501 antinucleon->setType(list[1]->getType());
1502
1503 ParticleList finallist;
1504
1505 finallist.push_back(nucleon);
1506 finallist.push_back(antinucleon);
1507
1508 if(list.size() > 2){
1509 for (G4int i = 2; i < (G4int)(list.size()); i++) {
1510 finallist.push_back(list[i]);
1511 }
1512 }
1513
1514 if(finallist.size()==2){
1515 G4double mn=nucleon->getMass();
1516 G4double my=antinucleon->getMass();
1517
1518 G4double ey=(sqrtS*sqrtS+my*my-mn*mn)/(2*sqrtS);
1519 G4double en=std::sqrt(ey*ey-my*my+mn*mn);
1520 nucleon->setEnergy(en);
1521 antinucleon->setEnergy(ey);
1522 G4double py=std::sqrt(ey*ey-my*my);
1523
1524 ThreeVector mom_antinucleon = Random::normVector(py);
1525 antinucleon->setMomentum(mom_antinucleon);
1526 nucleon->setMomentum(-mom_antinucleon);
1527 }
1528 else if(finallist.size() > 2){
1529 PhaseSpaceGenerator::generate(sqrtS, finallist);
1530 }
1531 else{
1532 INCL_ERROR("less than 2 mesons in NNbar annihilation!" << '\n');
1533 }
1534
1535
1536 for (G4int i = 0; i < 2; i++) {
1537 fs->addModifiedParticle(finallist[i]);
1538 }
1539 if(finallist.size()>2){
1540 for (G4int i = 2; i < (G4int)(list.size()); i++) {
1541 fs->addCreatedParticle(finallist[i]);
1542 }
1543 }
1544
1545
1546 //fs->addDestroyedParticle(nucleon);
1547 //fs->addDestroyedParticle(antinucleon);
1548
1549
1550 }
1551}
1552
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition G4Exception.hh:40
INCL_ERROR
#define INCL_ERROR(x)
Definition G4INCLLogger.hh:225
INCL_DEBUG
#define INCL_DEBUG(x)
Definition G4INCLLogger.hh:240
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4INCL::Config::getINCLXXDataFilePath
std::string const & getINCLXXDataFilePath() const
Set the ABLAXX datafile path.
Definition G4INCLConfig.hh:234
G4INCL::FinalState::addModifiedParticle
void addModifiedParticle(Particle *p)
Definition G4INCLFinalState.cc:60
G4INCL::FinalState::addCreatedParticle
void addCreatedParticle(Particle *p)
Definition G4INCLFinalState.cc:75
G4INCL::NNbarToAnnihilationChannel::findStringNumber
G4int findStringNumber(G4double rdm, std::vector< G4double > yields)
Definition G4INCLNNbarToAnnihilationChannel.cc:80
G4INCL::NNbarToAnnihilationChannel::read_file
G4double read_file(std::string filename, std::vector< G4double > &probabilities, std::vector< std::vector< std::string > > &particle_types)
Definition G4INCLNNbarToAnnihilationChannel.cc:56
G4INCL::NNbarToAnnihilationChannel::NNbarToAnnihilationChannel
NNbarToAnnihilationChannel(Nucleus *, Particle *, Particle *)
Definition G4INCLNNbarToAnnihilationChannel.cc:49
G4INCL::KinematicsUtils::compute_xs
G4double compute_xs(const std::vector< G4double > coefficients, const G4double pLab)
Definition G4INCLKinematicsUtils.cc:49
G4INCL::KinematicsUtils::totalEnergyInCM
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
Definition G4INCLKinematicsUtils.cc:121
G4INCL::KinematicsUtils::momentumInLab
G4double momentumInLab(Particle const *const p1, Particle const *const p2)
gives the momentum in the lab frame of two particles.
Definition G4INCLKinematicsUtils.cc:162
G4INCL::PhaseSpaceGenerator::generate
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
Definition G4INCLPhaseSpaceGenerator.cc:94
G4INCL::Random::normVector
ThreeVector normVector(G4double norm=1.)
Definition G4INCLRandom.cc:162
G4INCL::Random::shoot
G4double shoot()
Definition G4INCLRandom.cc:93