| 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 #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
47 namespace 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
56 G4double 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
80 G4int 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
102 void 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 G4double rdm = Random::shoot();
119
120 const std::vector<G4double> BFMM6 = {66.098, 0.153, -4.576, -38.319, 6.625}; //ppbar annihilation xs
121 const std::vector<G4double> BFMM1 = {119.066, 6.251, -0.006, -60.046, 11.958}; //ppbar total xs
122 const std::vector<G4double> BFMM471 = {108.104, 15.708, 0.832, -54.632, -6.958}; //npbar total xs
123
124 //PPbar annihilation xs
125 const std::vector<G4double> PPbar_pip_pim = {0.637, -0.340, -0.003, -0.439, 0.144};
126 const std::vector<G4double> PPbar_pip_pim_pi0 = {-2.065, 4.893, -1.130, 1.231, -0.212};
127 const std::vector<G4double> PPbar_pip_pim_omega = {3.020, 0.425, -0.029, -3.420, 0.867};
128 const std::vector<G4double> PPbar_pip_pim_Kp_Km = {-1.295, 1.897, -0.001, -0.365, 0.044};
129 const std::vector<G4double> PPbar_pip_pim_pi0_Kp_Km = {-12.220, 12.509, -0.351, 4.682, -0.777};
130 const std::vector<G4double> PPbar_2pip_2pim = {3.547, 0.095, 0.957, -3.444, 0.685};
131 const std::vector<G4double> PPbar_2pip_2pim_pi0 = {13.044, 1.449, 0.695, -12.313, 1.627};
132 const std::vector<G4double> PPbar_2pip_2pim_3pi0 = {6.398, 0.199, -1.103, -1.271, -0.380};
133 const std::vector<G4double> PPbar_3pip_3pim = {1.490, 0.240, 0.002, -1.012, 0.134};
134 const std::vector<G4double> PPbar_3pip_3pim_pi0 = {0.286, 1.634, -1.369, 3.099, -1.294};
135 const std::vector<G4double> PPbar_3pip_3pim_2pi0 = {-11.370, 12.503, -0.680, 10.059, -2.501};
136 const std::vector<G4double> PPbar_3pip_3pim_3pi0 = {-14.732, 12.338, -0.724, 11.342, -2.224};
137 const std::vector<G4double> PPbar_4pip_4pim = {-1.574, 1.607, -0.864, 1.253, -0.276};
138 const std::vector<G4double> PPbar_4pip_4pim_pi0 = {-1.096, 0.977, -0.995, 1.007, -0.171};
139
140 //NPbar annihilation xs
141 const std::vector<G4double> NPbar_pip_2pim = {-12.116, 14.485, -0.094, -1.632, 0.882, 5.000};
142 const std::vector<G4double> NPbar_pip_2pim_2pi0 = {8.276, 5.057, 0.483, -15.864, 2.552, 7.000};
143 const std::vector<G4double> NPbar_pip_2pim_3pi0 = {-1.500, 9.574, 0.528, -11.633, -0.615, 7.000};
144 const std::vector<G4double> NPbar_pip_2pim_pi0 = {7.999, 4.135, 0.608, -14.136, 1.590, 7.000};
145 const std::vector<G4double> NPbar_pip_pim_pi0_Km_K0 = {0.083, 0.091, -1.709, 0.284, -0.107};
146 const std::vector<G4double> NPbar_pip_pim_Km_K0 = {0.003, 0.297, -0.001, -0.143, 0.052};
147 const std::vector<G4double> NPbar_2pip_3pim_pi0 = {-14.701, 22.258, -0.001, -3.094, -0.190};
148 const std::vector<G4double> NPbar_2pip_3pim = {-0.616, 4.575, -0.002, -1.921, -0.153};
149
150
151 // File names
152 #ifdef INCLXX_IN_GEANT4_MODE
153 if(!G4FindDataDir("G4INCLDATA")) {
154 G4ExceptionDescription ed;
155 ed << " Data missing: set environment variable G4INCLDATA\n"
156 << " to point to the directory containing data files needed\n"
157 << " by the INCL++ model" << G4endl;
158 G4Exception("G4INCLDataFile::readData()","inflightppbarFS.dat, ...",
159 FatalException, ed);
160 }
161 G4String dataPath0{G4FindDataDir("G4INCLDATA")};
162 G4String dataPathppbar(dataPath0 + "/inflightppbarFS.dat");
163 G4String dataPathnpbar(dataPath0 + "/inflightnpbarFS.dat");
164 G4String dataPathppbark(dataPath0 + "/inflightppbarFSkaonic.dat");
165 G4String dataPathnpbark(dataPath0 + "/inflightnpbarFSkaonic.dat");
166 G4String dataPathpnbar(dataPath0 + "/inflightpnbarFS.dat"); //nbar case
167 G4String dataPathpnbark(dataPath0 + "/inflightpnbarFSkaonic.dat"); // nbar case
168 #else
169 //Config *theConfig = new G4INCL::Config;
170 //theConfig->setINCLXXDataFilePath(G4INCL::theINCLXXDataFilePath);
171 Config const *theConfig=theNucleus->getStore()->getConfig();
172 std::string path;
173 if(theConfig)
174 path = theConfig->getINCLXXDataFilePath();
175 std::string dataPathppbar(path + "/inflightppbarFS.dat");
176 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
177 std::string dataPathnpbar(path + "/inflightnpbarFS.dat");
178 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
179 std::string dataPathppbark(path + "/inflightppbarFSkaonic.dat");
180 INCL_DEBUG("Reading https://doi.org/10.1016/j.physrep.2005.03.002 ppbar kaonic final states" << dataPathppbark << '\n');
181 std::string dataPathnpbark(path + "/inflightnpbarFSkaonic.dat");
182 INCL_DEBUG("Reading https://doi.org/10.1007/BF02818764 and https://link.springer.com/article/10.1007/BF02754930 npbar kaonic final states" << dataPathnpbark << '\n');
183 std::string dataPathpnbar(path + "/inflightpnbarFS.dat"); // nbar case
184 std::string dataPathpnbark(path + "/inflightpnbarFSkaonic.dat"); // nbar case
185 #endif
186 /*std::string path = {"/home/zdemid/INCL/inclcode/data"};
187 std::string dataPathppbar(path + "/inflightppbarFS.dat");
188 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
189 std::string dataPathnpbar(path + "/inflightnpbarFS.dat");
190 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
191 std::string dataPathppbark(path + "/inflightppbarFSkaonic.dat");
192 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar kaonic final states" << dataPathppbark << '\n');
193 std::string dataPathnpbark(path + "/inflightnpbarFSkaonic.dat");
194 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar kaonic final states" << dataPathnpbark << '\n');
195 every time we remove lines for which we have data from BFMM
196 */
197
198 std::vector<G4double> probabilities; //will store each FS yield
199 std::vector<std::vector<std::string>> particle_types; //will store particle names
200 G4double sum; //will contain a sum of probabilities of all FS in the file
201 const G4double kaonicFSprob=0.05; //probability to kave kaonic FS
202
203
204 ParticleList list;
205 //list.push_back(nucleon);
206 //list.push_back(antinucleon);
207 // NNbar will not be in the list because they annihilate
208 const ThreeVector &rcol = nucleon->getPosition();
209 const ThreeVector zero;
210
211 //setting types of new particles and pushing them back to the list
212 if(nucleon->getType()==Neutron && antinucleon->getType()==antiProton){
213 //std::cout << "npbar"<< std::endl;
214 const G4double totalpnbar = KinematicsUtils::compute_xs(BFMM6, plab)*KinematicsUtils::compute_xs(BFMM471, plab)/KinematicsUtils::compute_xs(BFMM1, plab);
215 // xs is same for npbar, but the fs has different charge
216
217 if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab)) {
218 // First condition
219 Particle* p1 = new Particle(PiPlus, zero, rcol);
220 Particle* p2 = new Particle(PiMinus, zero, rcol);
221 Particle* p3 = new Particle(PiMinus, zero, rcol);
222
223 list.push_back(p1);
224 list.push_back(p2);
225 list.push_back(p3);
226 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
227 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab)) {
228 // Second condition
229 Particle* p1 = new Particle(PiPlus, zero, rcol);
230 Particle* p2 = new Particle(PiMinus, zero, rcol);
231 Particle* p3 = new Particle(PiZero, zero, rcol);
232 Particle* p4 = new Particle(PiZero, zero, rcol);
233 Particle* p5 = new Particle(PiMinus, zero, rcol);
234
235 list.push_back(p1);
236 list.push_back(p2);
237 list.push_back(p3);
238 list.push_back(p4);
239 list.push_back(p5);
240 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
241 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
242 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab)) {
243 // Third condition
244 Particle* p1 = new Particle(PiPlus, zero, rcol);
245 Particle* p2 = new Particle(PiMinus, zero, rcol);
246 Particle* p3 = new Particle(PiZero, zero, rcol);
247 Particle* p4 = new Particle(PiZero, zero, rcol);
248 Particle* p5 = new Particle(PiZero, zero, rcol);
249 Particle* p6 = new Particle(PiMinus, zero, rcol);
250
251 list.push_back(p1);
252 list.push_back(p2);
253 list.push_back(p3);
254 list.push_back(p4);
255 list.push_back(p5);
256 list.push_back(p6);
257 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
258 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
259 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
260 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab)) {
261 // Fourth condition
262 Particle* p1 = new Particle(PiPlus, zero, rcol);
263 Particle* p2 = new Particle(PiMinus, zero, rcol);
264 Particle* p3 = new Particle(PiZero, zero, rcol);
265 Particle* p4 = new Particle(PiMinus, zero, rcol);
266
267 list.push_back(p1);
268 list.push_back(p2);
269 list.push_back(p3);
270 list.push_back(p4);
271 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
272 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
273 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
274 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
275 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab)) {
276 // Fifth condition
277 Particle* p1 = new Particle(PiPlus, zero, rcol);
278 Particle* p2 = new Particle(PiMinus, zero, rcol);
279 Particle* p3 = new Particle(PiZero, zero, rcol);
280 Particle* p4 = new Particle(KMinus, zero, rcol);
281 Particle* p5 = new Particle(KZero, 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 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
289 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
290 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
291 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
292 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
293 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab)) {
294 // Sixth condition
295 Particle* p1 = new Particle(PiPlus, zero, rcol);
296 Particle* p2 = new Particle(PiMinus, zero, rcol);
297 Particle* p3 = new Particle(KMinus, zero, rcol);
298 Particle* p4 = new Particle(KZero, zero, rcol);
299
300 list.push_back(p1);
301 list.push_back(p2);
302 list.push_back(p3);
303 list.push_back(p4);
304 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
305 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
306 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
307 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
308 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
309 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
310 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab)) {
311 // Seventh condition
312 Particle* p1 = new Particle(PiPlus, zero, rcol);
313 Particle* p2 = new Particle(PiPlus, zero, rcol);
314 Particle* p3 = new Particle(PiMinus, zero, rcol);
315 Particle* p4 = new Particle(PiMinus, zero, rcol);
316 Particle* p5 = new Particle(PiMinus, zero, rcol);
317 Particle* p6 = new Particle(PiZero, zero, rcol);
318
319 list.push_back(p1);
320 list.push_back(p2);
321 list.push_back(p3);
322 list.push_back(p4);
323 list.push_back(p5);
324 list.push_back(p6);
325 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
326 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
327 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
328 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
329 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
330 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
331 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab) +
332 KinematicsUtils::compute_xs(NPbar_2pip_3pim, plab)) {
333 // Eighth condition
334 Particle* p1 = new Particle(PiPlus, zero, rcol);
335 Particle* p2 = new Particle(PiPlus, zero, rcol);
336 Particle* p3 = new Particle(PiMinus, zero, rcol);
337 Particle* p4 = new Particle(PiMinus, zero, rcol);
338 Particle* p5 = new Particle(PiMinus, zero, rcol);
339
340 list.push_back(p1);
341 list.push_back(p2);
342 list.push_back(p3);
343 list.push_back(p4);
344 list.push_back(p5);
345 } else {
346 // Default condition
347 //std::cout << "default condition pnbar"<< std::endl;
348 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
349 INCL_DEBUG("pionic npbar final state chosen" << '\n');
350 sum = read_file(dataPathnpbar, probabilities, particle_types);
351 rdm = (rdm/(1.-kaonicFSprob))*sum; //normalize by the sum of probabilities in the file
352 //now get the line number in the file where the FS particles are stored:
353 G4int n = findStringNumber(rdm, probabilities)-1;
354 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
355 if(particle_types[n][j] == "pi0"){
356 Particle *p = new Particle(PiZero, zero, rcol);
357 list.push_back(p);
358 }
359 else if(particle_types[n][j] == "pi-"){
360 Particle *p = new Particle(PiMinus, zero, rcol);
361 list.push_back(p);
362 }
363 else if(particle_types[n][j] == "pi+"){
364 Particle *p = new Particle(PiPlus, zero, rcol);
365 list.push_back(p);
366 }
367 else if(particle_types[n][j] == "omega"){
368 Particle *p = new Particle(Omega, zero, rcol);
369 list.push_back(p);
370 }
371 else if(particle_types[n][j] == "eta"){
372 Particle *p = new Particle(Eta, zero, rcol);
373 list.push_back(p);
374 }
375 else{
376 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
377 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
378 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
379 }
380 }
381 }
382 } // end of pionic option
383 else{
384 INCL_DEBUG("kaonic npbar final state chosen" << '\n');
385 sum = read_file(dataPathnpbark, probabilities, particle_types);
386 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
387 //now get the line number in the file where the FS particles are stored:
388 G4int n = findStringNumber(rdm, probabilities)-1;
389 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
390 if(particle_types[n][j] == "pi0"){
391 Particle *p = new Particle(PiZero, zero, rcol);
392 list.push_back(p);
393 }
394 else if(particle_types[n][j] == "pi-"){
395 Particle *p = new Particle(PiMinus, zero, rcol);
396 list.push_back(p);
397 }
398 else if(particle_types[n][j] == "pi+"){
399 Particle *p = new Particle(PiPlus, zero, rcol);
400 list.push_back(p);
401 }
402 else if(particle_types[n][j] == "omega"){
403 Particle *p = new Particle(Omega, zero, rcol);
404 list.push_back(p);
405 }
406 else if(particle_types[n][j] == "eta"){
407 Particle *p = new Particle(Eta, zero, rcol);
408 list.push_back(p);
409 }
410 else if(particle_types[n][j] == "K-"){
411 Particle *p = new Particle(KMinus, zero, rcol);
412 list.push_back(p);
413 }
414 else if(particle_types[n][j] == "K+"){
415 Particle *p = new Particle(KPlus, zero, rcol);
416 list.push_back(p);
417 }
418 else if(particle_types[n][j] == "K0"){
419 Particle *p = new Particle(KZero, zero, rcol);
420 list.push_back(p);
421 }
422 else if(particle_types[n][j] == "K0b"){
423 Particle *p = new Particle(KZeroBar, zero, rcol);
424 list.push_back(p);
425 }
426 else{
427 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
428 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
429 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
430 }
431 }
432 }
433 } // end of kaonic option
434 } // end of default annihilation
435
436 }
437 else if(nucleon->getType()==Proton && antinucleon->getType()==antiNeutron){
438 const G4double totalpnbar = KinematicsUtils::compute_xs(BFMM6, plab)*KinematicsUtils::compute_xs(BFMM471, plab)/KinematicsUtils::compute_xs(BFMM1, plab);
439 // xs is same for npbar, but the fs has different charge
440
441 if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab)) {
442 // First condition
443 Particle* p1 = new Particle(PiPlus, zero, rcol);
444 Particle* p2 = new Particle(PiMinus, zero, rcol);
445 Particle* p3 = new Particle(PiPlus, zero, rcol);
446
447 list.push_back(p1);
448 list.push_back(p2);
449 list.push_back(p3);
450 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
451 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab)) {
452 // Second condition
453 Particle* p1 = new Particle(PiPlus, zero, rcol);
454 Particle* p2 = new Particle(PiMinus, zero, rcol);
455 Particle* p3 = new Particle(PiZero, zero, rcol);
456 Particle* p4 = new Particle(PiZero, zero, rcol);
457 Particle* p5 = new Particle(PiPlus, zero, rcol);
458
459 list.push_back(p1);
460 list.push_back(p2);
461 list.push_back(p3);
462 list.push_back(p4);
463 list.push_back(p5);
464 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
465 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
466 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab)) {
467 // Third condition
468 Particle* p1 = new Particle(PiPlus, zero, rcol);
469 Particle* p2 = new Particle(PiMinus, zero, rcol);
470 Particle* p3 = new Particle(PiZero, zero, rcol);
471 Particle* p4 = new Particle(PiZero, zero, rcol);
472 Particle* p5 = new Particle(PiZero, zero, rcol);
473 Particle* p6 = new Particle(PiPlus, zero, rcol);
474
475 list.push_back(p1);
476 list.push_back(p2);
477 list.push_back(p3);
478 list.push_back(p4);
479 list.push_back(p5);
480 list.push_back(p6);
481 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
482 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
483 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
484 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab)) {
485 // Fourth condition
486 Particle* p1 = new Particle(PiPlus, zero, rcol);
487 Particle* p2 = new Particle(PiMinus, zero, rcol);
488 Particle* p3 = new Particle(PiZero, zero, rcol);
489 Particle* p4 = new Particle(PiPlus, zero, rcol);
490
491 list.push_back(p1);
492 list.push_back(p2);
493 list.push_back(p3);
494 list.push_back(p4);
495 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
496 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
497 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
498 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
499 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab)) {
500 // Fifth condition
501 Particle* p1 = new Particle(PiPlus, zero, rcol);
502 Particle* p2 = new Particle(PiMinus, zero, rcol);
503 Particle* p3 = new Particle(PiZero, zero, rcol);
504 Particle* p4 = new Particle(KPlus, zero, rcol);
505 Particle* p5 = new Particle(KZeroBar, zero, rcol);
506
507 list.push_back(p1);
508 list.push_back(p2);
509 list.push_back(p3);
510 list.push_back(p4);
511 list.push_back(p5);
512 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
513 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
514 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
515 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
516 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
517 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab)) {
518 // Sixth condition
519 Particle* p1 = new Particle(PiPlus, zero, rcol);
520 Particle* p2 = new Particle(PiMinus, zero, rcol);
521 Particle* p3 = new Particle(KPlus, zero, rcol);
522 Particle* p4 = new Particle(KZeroBar, zero, rcol);
523
524 list.push_back(p1);
525 list.push_back(p2);
526 list.push_back(p3);
527 list.push_back(p4);
528 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
529 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
530 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
531 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
532 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
533 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
534 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab)) {
535 // Seventh condition
536 Particle* p1 = new Particle(PiPlus, zero, rcol);
537 Particle* p2 = new Particle(PiPlus, zero, rcol);
538 Particle* p3 = new Particle(PiMinus, zero, rcol);
539 Particle* p4 = new Particle(PiMinus, zero, rcol);
540 Particle* p5 = new Particle(PiPlus, zero, rcol);
541 Particle* p6 = new Particle(PiZero, zero, rcol);
542
543 list.push_back(p1);
544 list.push_back(p2);
545 list.push_back(p3);
546 list.push_back(p4);
547 list.push_back(p5);
548 list.push_back(p6);
549 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
550 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
551 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
552 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
553 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
554 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
555 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab) +
556 KinematicsUtils::compute_xs(NPbar_2pip_3pim, plab)) {
557 // Eighth condition
558 Particle* p1 = new Particle(PiPlus, zero, rcol);
559 Particle* p2 = new Particle(PiPlus, zero, rcol);
560 Particle* p3 = new Particle(PiMinus, zero, rcol);
561 Particle* p4 = new Particle(PiMinus, zero, rcol);
562 Particle* p5 = new Particle(PiPlus, zero, rcol);
563
564 list.push_back(p1);
565 list.push_back(p2);
566 list.push_back(p3);
567 list.push_back(p4);
568 list.push_back(p5);
569 } else {
570 // Default condition
571 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
572 INCL_DEBUG("pionic pnbar final state chosen" << '\n');
573 sum = read_file(dataPathpnbar, probabilities, particle_types);
574 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.95 normalize by the sum of probabilities in the file
575 //now get the line number in the file where the FS particles are stored:
576 G4int n = findStringNumber(rdm, probabilities)-1;
577 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
578 if(particle_types[n][j] == "pi0"){
579 Particle *p = new Particle(PiZero, zero, rcol);
580 list.push_back(p);
581 }
582 else if(particle_types[n][j] == "pi-"){
583 Particle *p = new Particle(PiMinus, zero, rcol);
584 list.push_back(p);
585 }
586 else if(particle_types[n][j] == "pi+"){
587 Particle *p = new Particle(PiPlus, zero, rcol);
588 list.push_back(p);
589 }
590 else if(particle_types[n][j] == "omega"){
591 Particle *p = new Particle(Omega, zero, rcol);
592 list.push_back(p);
593 }
594 else if(particle_types[n][j] == "eta"){
595 Particle *p = new Particle(Eta, zero, rcol);
596 list.push_back(p);
597 }
598 else{
599 INCL_ERROR("Some non-existing FS particle detected when reading nbar FS files");
600 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
601 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
602 }
603 }
604 }
605 } // end of pionic option
606 else{
607 INCL_DEBUG("kaonic pnbar final state chosen" << '\n');
608 sum = read_file(dataPathnpbark, probabilities, particle_types);
609 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
610 //now get the line number in the file where the FS particles are stored:
611 G4int n = findStringNumber(rdm, probabilities)-1;
612 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
613 if(particle_types[n][j] == "pi0"){
614 Particle *p = new Particle(PiZero, zero, rcol);
615 list.push_back(p);
616 }
617 else if(particle_types[n][j] == "pi-"){
618 Particle *p = new Particle(PiMinus, zero, rcol);
619 list.push_back(p);
620 }
621 else if(particle_types[n][j] == "pi+"){
622 Particle *p = new Particle(PiPlus, zero, rcol);
623 list.push_back(p);
624 }
625 else if(particle_types[n][j] == "omega"){
626 Particle *p = new Particle(Omega, zero, rcol);
627 list.push_back(p);
628 }
629 else if(particle_types[n][j] == "eta"){
630 Particle *p = new Particle(Eta, zero, rcol);
631 list.push_back(p);
632 }
633 else if(particle_types[n][j] == "K-"){
634 Particle *p = new Particle(KMinus, zero, rcol);
635 list.push_back(p);
636 }
637 else if(particle_types[n][j] == "K+"){
638 Particle *p = new Particle(KPlus, zero, rcol);
639 list.push_back(p);
640 }
641 else if(particle_types[n][j] == "K0"){
642 Particle *p = new Particle(KZero, zero, rcol);
643 list.push_back(p);
644 }
645 else if(particle_types[n][j] == "K0b"){
646 Particle *p = new Particle(KZeroBar, zero, rcol);
647 list.push_back(p);
648 }
649 else{
650 INCL_ERROR("Some non-existing FS particle detected when reading pnbar FS files");
651 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
652 std::cout << "gotcha! " << particle_types[n][jj] << std::endl; }
653 }
654 }
655 } // end of kaonic option
656 } // end of default annihilation
657
658 }
659 else{ //ppbar or nnbar
660 //std::cout << "ppbar or nnbar"<< std::endl;
661 const G4double totalppbar = KinematicsUtils::compute_xs(BFMM6, plab);
662 // same for nnbar
663
664 if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab)) {
665 // First condition
666 Particle* p1 = new Particle(PiPlus, zero, rcol);
667 Particle* p2 = new Particle(PiMinus, zero, rcol);
668
669 list.push_back(p1);
670 list.push_back(p2);
671
672
673 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
674 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab)) {
675 // Second condition
676 Particle* p1 = new Particle(PiPlus, zero, rcol);
677 Particle* p2 = new Particle(PiMinus, zero, rcol);
678 Particle* p3 = new Particle(PiZero, zero, rcol);
679
680 list.push_back(p1);
681 list.push_back(p2);
682 list.push_back(p3);
683 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
684 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
685 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab)) {
686 // Third condition
687 Particle* p1 = new Particle(PiPlus, zero, rcol);
688 Particle* p2 = new Particle(PiMinus, zero, rcol);
689 Particle* p3 = new Particle(Omega, zero, rcol);
690
691 list.push_back(p1);
692 list.push_back(p2);
693 list.push_back(p3);
694 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
695 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
696 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
697 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab)) {
698 // Fourth condition
699 Particle* p1 = new Particle(PiPlus, zero, rcol);
700 Particle* p2 = new Particle(PiMinus, zero, rcol);
701 Particle* p3 = new Particle(KPlus, zero, rcol);
702 Particle* p4 = new Particle(KMinus, zero, rcol);
703
704 list.push_back(p1);
705 list.push_back(p2);
706 list.push_back(p3);
707 list.push_back(p4);
708 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
709 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
710 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
711 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
712 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab)) {
713 // Fifth condition
714 Particle* p1 = new Particle(PiPlus, zero, rcol);
715 Particle* p2 = new Particle(PiMinus, zero, rcol);
716 Particle* p3 = new Particle(PiZero, zero, rcol);
717 Particle* p4 = new Particle(KPlus, zero, rcol);
718 Particle* p5 = new Particle(KMinus, zero, rcol);
719
720 list.push_back(p1);
721 list.push_back(p2);
722 list.push_back(p3);
723 list.push_back(p4);
724 list.push_back(p5);
725 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
726 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
727 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
728 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
729 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
730 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab)) {
731 // Sixth condition
732 Particle* p1 = new Particle(PiPlus, zero, rcol);
733 Particle* p2 = new Particle(PiMinus, zero, rcol);
734 Particle* p3 = new Particle(PiPlus, zero, rcol);
735 Particle* p4 = new Particle(PiMinus, zero, rcol);
736
737 list.push_back(p1);
738 list.push_back(p2);
739 list.push_back(p3);
740 list.push_back(p4);
741 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
742 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
743 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
744 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
745 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
746 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
747 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab)) {
748 // Seventh condition
749 Particle* p1 = new Particle(PiPlus, zero, rcol);
750 Particle* p2 = new Particle(PiMinus, zero, rcol);
751 Particle* p3 = new Particle(PiPlus, zero, rcol);
752 Particle* p4 = new Particle(PiMinus, zero, rcol);
753 Particle* p5 = new Particle(PiZero, zero, rcol);
754
755 list.push_back(p1);
756 list.push_back(p2);
757 list.push_back(p3);
758 list.push_back(p4);
759 list.push_back(p5);
760 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
761 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
762 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
763 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
764 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
765 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
766 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
767 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab)) {
768 // Eighth condition
769 Particle* p1 = new Particle(PiPlus, zero, rcol);
770 Particle* p2 = new Particle(PiMinus, zero, rcol);
771 Particle* p3 = new Particle(PiPlus, zero, rcol);
772 Particle* p4 = new Particle(PiMinus, zero, rcol);
773 Particle* p5 = new Particle(PiZero, zero, rcol);
774 Particle* p6 = new Particle(PiZero, zero, rcol);
775 Particle* p7 = new Particle(PiZero, zero, rcol);
776
777 list.push_back(p1);
778 list.push_back(p2);
779 list.push_back(p3);
780 list.push_back(p4);
781 list.push_back(p5);
782 list.push_back(p6);
783 list.push_back(p7);
784 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
785 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
786 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
787 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
788 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
789 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
790 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
791 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
792 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab)) {
793 // Ninth condition
794 Particle* p1 = new Particle(PiPlus, zero, rcol);
795 Particle* p2 = new Particle(PiMinus, zero, rcol);
796 Particle* p3 = new Particle(PiPlus, zero, rcol);
797 Particle* p4 = new Particle(PiMinus, zero, rcol);
798 Particle* p5 = new Particle(PiPlus, zero, rcol);
799 Particle* p6 = new Particle(PiMinus, zero, rcol);
800
801 list.push_back(p1);
802 list.push_back(p2);
803 list.push_back(p3);
804 list.push_back(p4);
805 list.push_back(p5);
806 list.push_back(p6);
807 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
808 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
809 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
810 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
811 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
812 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
813 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
814 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
815 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
816 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab)) {
817 // Tenth condition
818 Particle* p1 = new Particle(PiPlus, zero, rcol);
819 Particle* p2 = new Particle(PiMinus, zero, rcol);
820 Particle* p3 = new Particle(PiPlus, zero, rcol);
821 Particle* p4 = new Particle(PiMinus, zero, rcol);
822 Particle* p5 = new Particle(PiPlus, zero, rcol);
823 Particle* p6 = new Particle(PiMinus, zero, rcol);
824 Particle* p7 = new Particle(PiZero, zero, rcol);
825
826 list.push_back(p1);
827 list.push_back(p2);
828 list.push_back(p3);
829 list.push_back(p4);
830 list.push_back(p5);
831 list.push_back(p6);
832 list.push_back(p7);
833 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
834 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
835 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
836 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
837 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
838 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
839 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
840 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
841 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
842 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
843 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab)) {
844 // Eleventh condition
845 Particle* p1 = new Particle(PiPlus, zero, rcol);
846 Particle* p2 = new Particle(PiMinus, zero, rcol);
847 Particle* p3 = new Particle(PiPlus, zero, rcol);
848 Particle* p4 = new Particle(PiMinus, zero, rcol);
849 Particle* p5 = new Particle(PiPlus, zero, rcol);
850 Particle* p6 = new Particle(PiMinus, zero, rcol);
851 Particle* p7 = new Particle(PiZero, zero, rcol);
852 Particle* p8 = new Particle(PiZero, zero, rcol);
853
854 list.push_back(p1);
855 list.push_back(p2);
856 list.push_back(p3);
857 list.push_back(p4);
858 list.push_back(p5);
859 list.push_back(p6);
860 list.push_back(p7);
861 list.push_back(p8);
862 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
863 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
864 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
865 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
866 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
867 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
868 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
869 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
870 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
871 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
872 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
873 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab)) {
874 // Twelfth condition
875 Particle* p1 = new Particle(PiPlus, zero, rcol);
876 Particle* p2 = new Particle(PiMinus, zero, rcol);
877 Particle* p3 = new Particle(PiPlus, zero, rcol);
878 Particle* p4 = new Particle(PiMinus, zero, rcol);
879 Particle* p5 = new Particle(PiPlus, zero, rcol);
880 Particle* p6 = new Particle(PiMinus, zero, rcol);
881 Particle* p7 = new Particle(PiZero, zero, rcol);
882 Particle* p8 = new Particle(PiZero, zero, rcol);
883 Particle* p9 = new Particle(PiZero, zero, rcol);
884
885 list.push_back(p1);
886 list.push_back(p2);
887 list.push_back(p3);
888 list.push_back(p4);
889 list.push_back(p5);
890 list.push_back(p6);
891 list.push_back(p7);
892 list.push_back(p8);
893 list.push_back(p9);
894 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
895 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
896 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
897 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
898 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
899 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
900 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
901 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
902 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
903 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
904 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
905 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab) +
906 KinematicsUtils::compute_xs(PPbar_4pip_4pim, plab)) {
907 // Thirteenth condition
908 Particle* p1 = new Particle(PiPlus, zero, rcol);
909 Particle* p2 = new Particle(PiMinus, zero, rcol);
910 Particle* p3 = new Particle(PiPlus, zero, rcol);
911 Particle* p4 = new Particle(PiMinus, zero, rcol);
912 Particle* p5 = new Particle(PiPlus, zero, rcol);
913 Particle* p6 = new Particle(PiMinus, zero, rcol);
914 Particle* p7 = new Particle(PiPlus, zero, rcol);
915 Particle* p8 = new Particle(PiMinus, zero, rcol);
916
917 list.push_back(p1);
918 list.push_back(p2);
919 list.push_back(p3);
920 list.push_back(p4);
921 list.push_back(p5);
922 list.push_back(p6);
923 list.push_back(p7);
924 list.push_back(p8);
925 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
926 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
927 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
928 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
929 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
930 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
931 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
932 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
933 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
934 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
935 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
936 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab) +
937 KinematicsUtils::compute_xs(PPbar_4pip_4pim, plab) +
938 KinematicsUtils::compute_xs(PPbar_4pip_4pim_pi0, plab)) {
939 // Fourteenth condition
940 Particle* p1 = new Particle(PiPlus, zero, rcol);
941 Particle* p2 = new Particle(PiMinus, zero, rcol);
942 Particle* p3 = new Particle(PiPlus, zero, rcol);
943 Particle* p4 = new Particle(PiMinus, zero, rcol);
944 Particle* p5 = new Particle(PiPlus, zero, rcol);
945 Particle* p6 = new Particle(PiMinus, zero, rcol);
946 Particle* p7 = new Particle(PiPlus, zero, rcol);
947 Particle* p8 = new Particle(PiMinus, zero, rcol);
948 Particle* p9 = new Particle(PiZero, zero, rcol);
949
950 list.push_back(p1);
951 list.push_back(p2);
952 list.push_back(p3);
953 list.push_back(p4);
954 list.push_back(p5);
955 list.push_back(p6);
956 list.push_back(p7);
957 list.push_back(p8);
958 list.push_back(p9);
959 } else {
960 // Default condition
961 if(rdm < (1.-kaonicFSprob)){ // pionic FS was chosen
962 INCL_DEBUG("pionic pp final state chosen" << '\n');
963 sum = read_file(dataPathppbar, probabilities, particle_types);
964 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.88 normalize by the sum of probabilities in the file
965 //now get the line number in the file where the FS particles are stored:
966 G4int n = findStringNumber(rdm, probabilities)-1;
967 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
968 if(particle_types[n][j] == "pi0"){
969 Particle *p = new Particle(PiZero, zero, rcol);
970 list.push_back(p);
971 }
972 else if(particle_types[n][j] == "pi-"){
973 Particle *p = new Particle(PiMinus, zero, rcol);
974 list.push_back(p);
975 }
976 else if(particle_types[n][j] == "pi+"){
977 Particle *p = new Particle(PiPlus, zero, rcol);
978 list.push_back(p);
979 }
980 else if(particle_types[n][j] == "omega"){
981 Particle *p = new Particle(Omega, zero, rcol);
982 list.push_back(p);
983 }
984 else if(particle_types[n][j] == "eta"){
985 Particle *p = new Particle(Eta, zero, rcol);
986 list.push_back(p);
987 }
988 else{
989 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
990 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
991 std::cout << "gotcha! " << particle_types[n][jj] << std::endl; }
992 }
993 }
994 } //end of pionic option
995 else{
996 INCL_DEBUG("kaonic pp final state chosen" << '\n');
997 sum = read_file(dataPathppbark, probabilities, particle_types);
998 rdm = ((1-rdm)/kaonicFSprob)*sum;//2670 normalize by the sum of probabilities in the file
999 //now get the line number in the file where the FS particles are stored:
1000 G4int n = findStringNumber(rdm, probabilities)-1;
1001 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
1002 if(particle_types[n][j] == "pi0"){
1003 Particle *p = new Particle(PiZero, zero, rcol);
1004 list.push_back(p);
1005 }
1006 else if(particle_types[n][j] == "pi-"){
1007 Particle *p = new Particle(PiMinus, zero, rcol);
1008 list.push_back(p);
1009 }
1010 else if(particle_types[n][j] == "pi+"){
1011 Particle *p = new Particle(PiPlus, zero, rcol);
1012 list.push_back(p);
1013 }
1014 else if(particle_types[n][j] == "omega"){
1015 Particle *p = new Particle(Omega, zero, rcol);
1016 list.push_back(p);
1017 }
1018 else if(particle_types[n][j] == "eta"){
1019 Particle *p = new Particle(Eta, zero, rcol);
1020 list.push_back(p);
1021 }
1022 else if(particle_types[n][j] == "K-"){
1023 Particle *p = new Particle(KMinus, zero, rcol);
1024 list.push_back(p);
1025 }
1026 else if(particle_types[n][j] == "K+"){
1027 Particle *p = new Particle(KPlus, zero, rcol);
1028 list.push_back(p);
1029 }
1030 else if(particle_types[n][j] == "K0"){
1031 Particle *p = new Particle(KZero, zero, rcol);
1032 list.push_back(p);
1033 }
1034 else if(particle_types[n][j] == "K0b"){
1035 Particle *p = new Particle(KZeroBar, zero, rcol);
1036 list.push_back(p);
1037 }
1038 else{
1039 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
1040 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
1041 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
1042 }
1043 }
1044 }
1045 } // end of kaonic option
1046 } // end of default condition
1047 } // end of ppbar and nnbar case
1048
1049
1050
1051 nucleon->setType(list[0]->getType());
1052 antinucleon->setType(list[1]->getType());
1053
1054 ParticleList finallist;
1055
1056 finallist.push_back(nucleon);
1057 finallist.push_back(antinucleon);
1058
1059 if(list.size() > 2){
1060 for (G4int i = 2; i < (G4int)(list.size()); i++) {
1061 finallist.push_back(list[i]);
1062 }
1063 }
1064
1065 if(finallist.size()==2){
1066 G4double mn=nucleon->getMass();
1067 G4double my=antinucleon->getMass();
1068
1069 G4double ey=(sqrtS*sqrtS+my*my-mn*mn)/(2*sqrtS);
1070 G4double en=std::sqrt(ey*ey-my*my+mn*mn);
1071 nucleon->setEnergy(en);
1072 antinucleon->setEnergy(ey);
1073 G4double py=std::sqrt(ey*ey-my*my);
1074
1075 ThreeVector mom_antinucleon = Random::normVector(py);
1076 antinucleon->setMomentum(mom_antinucleon);
1077 nucleon->setMomentum(-mom_antinucleon);
1078 }
1079 else if(finallist.size() > 2){
1080 PhaseSpaceGenerator::generate(sqrtS, finallist);
1081 }
1082 else{
1083 INCL_ERROR("less than 2 mesons in NNbar annihilation!" << '\n');
1084 }
1085
1086
1087 for (G4int i = 0; i < 2; i++) {
1088 fs->addModifiedParticle(finallist[i]);
1089 }
1090 if(finallist.size()>2){
1091 for (G4int i = 2; i < (G4int)(list.size()); i++) {
1092 fs->addCreatedParticle(finallist[i]);
1093 }
1094 }
1095
1096
1097 //fs->addDestroyedParticle(nucleon);
1098 //fs->addDestroyedParticle(antinucleon);
1099
1100
1101 }
1102 }
1103
1104