| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 #include "G4EnvironmentUtils.hh" 37 #include "G4EnvironmentUtils.hh"
38 #include "G4INCLNNbarToAnnihilationChannel.hh" 38 #include "G4INCLNNbarToAnnihilationChannel.hh"
39 #include "G4INCLKinematicsUtils.hh" 39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh" 40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh" 41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh" 42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh" 43 #include "G4INCLLogger.hh"
44 #include <algorithm> 44 #include <algorithm>
45 #include "G4INCLPhaseSpaceGenerator.hh" 45 #include "G4INCLPhaseSpaceGenerator.hh"
46 46
47 namespace G4INCL { 47 namespace G4INCL {
48 48
49 NNbarToAnnihilationChannel::NNbarToAnnihilat 49 NNbarToAnnihilationChannel::NNbarToAnnihilationChannel(Nucleus *n, Particle *p1, Particle *p2)
50 :theNucleus(n), particle1(p1), particle2(p 50 :theNucleus(n), particle1(p1), particle2(p2)
51 {} 51 {}
52 52
53 NNbarToAnnihilationChannel::~NNbarToAnnihila 53 NNbarToAnnihilationChannel::~NNbarToAnnihilationChannel(){}
54 54
55 //fill probabilities and particle types from d 55 //fill probabilities and particle types from datafile and return probability sum for normalization
56 G4double NNbarToAnnihilationChannel::read_file 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); 57 std::ifstream file(filename);
58 G4double sum_probs = 0.0; 58 G4double sum_probs = 0.0;
59 if (file.is_open()) { 59 if (file.is_open()) {
60 std::string line; 60 std::string line;
61 while (getline(file, line)) { 61 while (getline(file, line)) {
62 std::istringstream iss(line); 62 std::istringstream iss(line);
63 G4double prob; 63 G4double prob;
64 iss >> prob; 64 iss >> prob;
65 sum_probs += prob; 65 sum_probs += prob;
66 probabilities.push_back(prob); 66 probabilities.push_back(prob);
67 std::vector<std::string> types; 67 std::vector<std::string> types;
68 std::string type; 68 std::string type;
69 while (iss >> type) { 69 while (iss >> type) {
70 types.push_back(type); 70 types.push_back(type);
71 } 71 }
72 particle_types.push_back(std::move(t << 72 particle_types.push_back(types);
73 } 73 }
74 } 74 }
75 else std::cout << "ERROR no fread_file " << 75 else std::cout << "ERROR no fread_file " << filename << std::endl;
76 return sum_probs; 76 return sum_probs;
77 } 77 }
78 78
79 //this function will tell you the FS line numb 79 //this function will tell you the FS line number from the datafile based on your random value
80 G4int NNbarToAnnihilationChannel::findStringNu 80 G4int NNbarToAnnihilationChannel::findStringNumber(G4double rdm, std::vector<G4double> yields) {
81 G4int stringNumber = -1; 81 G4int stringNumber = -1;
82 G4double smallestsum = 0.0; 82 G4double smallestsum = 0.0;
83 G4double biggestsum = yields[0]; 83 G4double biggestsum = yields[0];
84 //std::cout << "initial input " << rdm << 84 //std::cout << "initial input " << rdm << std::endl;
85 for (G4int i = 0; i < static_cast<G4int>(y 85 for (G4int i = 0; i < static_cast<G4int>(yields.size()-1); i++) {
86 if (rdm >= smallestsum && rdm <= bigge 86 if (rdm >= smallestsum && rdm <= biggestsum) {
87 //std::cout << smallestsum << " an 87 //std::cout << smallestsum << " and " << biggestsum << std::endl;
88 stringNumber = i+1; 88 stringNumber = i+1;
89 } 89 }
90 smallestsum += yields[i]; 90 smallestsum += yields[i];
91 biggestsum += yields[i+1]; 91 biggestsum += yields[i+1];
92 } 92 }
93 if(stringNumber==-1) stringNumber = static 93 if(stringNumber==-1) stringNumber = static_cast<G4int>(yields.size());
94 if(stringNumber==-1){ 94 if(stringNumber==-1){
95 INCL_ERROR("ERROR in findStringNumber (s 95 INCL_ERROR("ERROR in findStringNumber (stringNumber=-1)");
96 std::cout << "ERROR in findStringNumber" 96 std::cout << "ERROR in findStringNumber" << std::endl;
97 } 97 }
98 return stringNumber; 98 return stringNumber;
99 } 99 }
100 100
101 101
102 void NNbarToAnnihilationChannel::fillFinalStat 102 void NNbarToAnnihilationChannel::fillFinalState(FinalState *fs) {
103 103
104 Particle *nucleon; 104 Particle *nucleon;
105 Particle *antinucleon; 105 Particle *antinucleon;
106 106
107 if(particle1->isNucleon()){ 107 if(particle1->isNucleon()){
108 nucleon = particle1; 108 nucleon = particle1;
109 antinucleon = particle2; 109 antinucleon = particle2;
110 } 110 }
111 else{ 111 else{
112 nucleon = particle2; 112 nucleon = particle2;
113 antinucleon = particle1; 113 antinucleon = particle1;
114 } 114 }
115 115
116 const G4double plab = 0.001*KinematicsUtil 116 const G4double plab = 0.001*KinematicsUtils::momentumInLab(particle1, particle2); //GeV
117 const G4double sqrtS = KinematicsUtils::to 117 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(nucleon, antinucleon);
118 G4double rdm = Random::shoot(); 118 G4double rdm = Random::shoot();
119 119
120 const std::vector<G4double> BFMM6 = {66.09 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.0 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 122 const std::vector<G4double> BFMM471 = {108.104, 15.708, 0.832, -54.632, -6.958}; //npbar total xs
123 123
124 //PPbar annihilation xs 124 //PPbar annihilation xs
125 const std::vector<G4double> PPbar_pip_pim 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_ 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_ 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_ 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_ 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_2pi 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_2pi 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_2pi 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_3pi 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_3pi 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_3pi 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_3pi 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_4pi 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_4pi 138 const std::vector<G4double> PPbar_4pip_4pim_pi0 = {-1.096, 0.977, -0.995, 1.007, -0.171};
139 139
140 //NPbar annihilation xs 140 //NPbar annihilation xs
141 const std::vector<G4double> NPbar_pip_2pim 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 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 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 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_ 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_ 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_3pi 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_3pi 148 const std::vector<G4double> NPbar_2pip_3pim = {-0.616, 4.575, -0.002, -1.921, -0.153};
149 149
150 150
151 // File names 151 // File names
152 #ifdef INCLXX_IN_GEANT4_MODE 152 #ifdef INCLXX_IN_GEANT4_MODE
153 if(!G4FindDataDir("G4INCLDATA")) { 153 if(!G4FindDataDir("G4INCLDATA")) {
154 G4ExceptionDescription ed; 154 G4ExceptionDescription ed;
155 ed << " Data missing: set environm 155 ed << " Data missing: set environment variable G4INCLDATA\n"
156 << " to point to the directory 156 << " to point to the directory containing data files needed\n"
157 << " by the INCL++ model" << G4 157 << " by the INCL++ model" << G4endl;
158 G4Exception("G4INCLDataFile::re 158 G4Exception("G4INCLDataFile::readData()","inflightppbarFS.dat, ...",
159 FatalException, ed); 159 FatalException, ed);
160 } 160 }
161 G4String dataPath0{G4FindDataDir("G4 161 G4String dataPath0{G4FindDataDir("G4INCLDATA")};
162 G4String dataPathppbar(dataPath0 + " 162 G4String dataPathppbar(dataPath0 + "/inflightppbarFS.dat");
163 G4String dataPathnpbar(dataPath0 + " 163 G4String dataPathnpbar(dataPath0 + "/inflightnpbarFS.dat");
164 G4String dataPathppbark(dataPath0 + 164 G4String dataPathppbark(dataPath0 + "/inflightppbarFSkaonic.dat");
165 G4String dataPathnpbark(dataPath0 + 165 G4String dataPathnpbark(dataPath0 + "/inflightnpbarFSkaonic.dat");
166 G4String dataPathpnbar(dataPath0 + " 166 G4String dataPathpnbar(dataPath0 + "/inflightpnbarFS.dat"); //nbar case
167 G4String dataPathpnbark(dataPath0 + 167 G4String dataPathpnbark(dataPath0 + "/inflightpnbarFSkaonic.dat"); // nbar case
168 #else 168 #else
169 //Config *theConfig = new G4INCL::Co 169 //Config *theConfig = new G4INCL::Config;
170 //theConfig->setINCLXXDataFilePath(G4I 170 //theConfig->setINCLXXDataFilePath(G4INCL::theINCLXXDataFilePath);
171 Config const *theConfig=theNucleus-> 171 Config const *theConfig=theNucleus->getStore()->getConfig();
172 std::string path; 172 std::string path;
173 if(theConfig) 173 if(theConfig)
174 path = theConfig->getINCLXXDataFil 174 path = theConfig->getINCLXXDataFilePath();
175 std::string dataPathppbar(path + "/i 175 std::string dataPathppbar(path + "/inflightppbarFS.dat");
176 INCL_DEBUG("Reading https://doi.org/ 176 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
177 std::string dataPathnpbar(path + "/i 177 std::string dataPathnpbar(path + "/inflightnpbarFS.dat");
178 INCL_DEBUG("Reading https://doi.org/ 178 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
179 std::string dataPathppbark(path + "/ 179 std::string dataPathppbark(path + "/inflightppbarFSkaonic.dat");
180 INCL_DEBUG("Reading https://doi.org/ 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 + "/ 181 std::string dataPathnpbark(path + "/inflightnpbarFSkaonic.dat");
182 INCL_DEBUG("Reading https://doi.org/ 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 + "/i 183 std::string dataPathpnbar(path + "/inflightpnbarFS.dat"); // nbar case
184 std::string dataPathpnbark(path + "/ 184 std::string dataPathpnbark(path + "/inflightpnbarFSkaonic.dat"); // nbar case
185 #endif 185 #endif
186 /*std::string path = {"/home/zdemid/IN 186 /*std::string path = {"/home/zdemid/INCL/inclcode/data"};
187 std::string dataPathppbar(path + "/inf 187 std::string dataPathppbar(path + "/inflightppbarFS.dat");
188 INCL_DEBUG("Reading https://doi.org/10 188 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
189 std::string dataPathnpbar(path + "/inf 189 std::string dataPathnpbar(path + "/inflightnpbarFS.dat");
190 INCL_DEBUG("Reading https://doi.org/10 190 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
191 std::string dataPathppbark(path + "/in 191 std::string dataPathppbark(path + "/inflightppbarFSkaonic.dat");
192 INCL_DEBUG("Reading https://doi.org/10 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 + "/in 193 std::string dataPathnpbark(path + "/inflightnpbarFSkaonic.dat");
194 INCL_DEBUG("Reading https://doi.org/10 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 w 195 every time we remove lines for which we have data from BFMM
196 */ 196 */
197 197
198 std::vector<G4double> probabilities; //wil 198 std::vector<G4double> probabilities; //will store each FS yield
199 std::vector<std::vector<std::string>> part 199 std::vector<std::vector<std::string>> particle_types; //will store particle names
200 G4double sum; //will contain a sum of prob 200 G4double sum; //will contain a sum of probabilities of all FS in the file
201 const G4double kaonicFSprob=0.05; //probab 201 const G4double kaonicFSprob=0.05; //probability to kave kaonic FS
202 202
203 203
204 ParticleList list; 204 ParticleList list;
205 //list.push_back(nucleon); 205 //list.push_back(nucleon);
206 //list.push_back(antinucleon); 206 //list.push_back(antinucleon);
207 // NNbar will not be in the list because t 207 // NNbar will not be in the list because they annihilate
208 const ThreeVector &rcol = nucleon->getPosi 208 const ThreeVector &rcol = nucleon->getPosition();
209 const ThreeVector zero; 209 const ThreeVector zero;
210 210
211 //setting types of new particles and pushi 211 //setting types of new particles and pushing them back to the list
212 if(nucleon->getType()==Neutron && antinucl 212 if(nucleon->getType()==Neutron && antinucleon->getType()==antiProton){
213 //std::cout << "npbar"<< std::endl; 213 //std::cout << "npbar"<< std::endl;
214 const G4double totalpnbar = KinematicsUt 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 215 // xs is same for npbar, but the fs has different charge
216 216
217 if (rdm * totalpnbar < KinematicsUtils:: 217 if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab)) {
218 // First condition 218 // First condition
219 Particle* p1 = new Particle(PiPlus, 219 Particle* p1 = new Particle(PiPlus, zero, rcol);
220 Particle* p2 = new Particle(PiMinus, 220 Particle* p2 = new Particle(PiMinus, zero, rcol);
221 Particle* p3 = new Particle(PiMinus, 221 Particle* p3 = new Particle(PiMinus, zero, rcol);
222 222
223 list.push_back(p1); 223 list.push_back(p1);
224 list.push_back(p2); 224 list.push_back(p2);
225 list.push_back(p3); 225 list.push_back(p3);
226 } else if (rdm * totalpnbar < Kinematics 226 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
227 KinematicsUt 227 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab)) {
228 // Second condition 228 // Second condition
229 Particle* p1 = new Particle(PiPlus, 229 Particle* p1 = new Particle(PiPlus, zero, rcol);
230 Particle* p2 = new Particle(PiMinus, 230 Particle* p2 = new Particle(PiMinus, zero, rcol);
231 Particle* p3 = new Particle(PiZero, 231 Particle* p3 = new Particle(PiZero, zero, rcol);
232 Particle* p4 = new Particle(PiZero, 232 Particle* p4 = new Particle(PiZero, zero, rcol);
233 Particle* p5 = new Particle(PiMinus, 233 Particle* p5 = new Particle(PiMinus, zero, rcol);
234 234
235 list.push_back(p1); 235 list.push_back(p1);
236 list.push_back(p2); 236 list.push_back(p2);
237 list.push_back(p3); 237 list.push_back(p3);
238 list.push_back(p4); 238 list.push_back(p4);
239 list.push_back(p5); 239 list.push_back(p5);
240 } else if (rdm * totalpnbar < Kinematics 240 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
241 KinematicsUt 241 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
242 KinematicsUt 242 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab)) {
243 // Third condition 243 // Third condition
244 Particle* p1 = new Particle(PiPlus, 244 Particle* p1 = new Particle(PiPlus, zero, rcol);
245 Particle* p2 = new Particle(PiMinus, 245 Particle* p2 = new Particle(PiMinus, zero, rcol);
246 Particle* p3 = new Particle(PiZero, 246 Particle* p3 = new Particle(PiZero, zero, rcol);
247 Particle* p4 = new Particle(PiZero, 247 Particle* p4 = new Particle(PiZero, zero, rcol);
248 Particle* p5 = new Particle(PiZero, 248 Particle* p5 = new Particle(PiZero, zero, rcol);
249 Particle* p6 = new Particle(PiMinus, 249 Particle* p6 = new Particle(PiMinus, zero, rcol);
250 250
251 list.push_back(p1); 251 list.push_back(p1);
252 list.push_back(p2); 252 list.push_back(p2);
253 list.push_back(p3); 253 list.push_back(p3);
254 list.push_back(p4); 254 list.push_back(p4);
255 list.push_back(p5); 255 list.push_back(p5);
256 list.push_back(p6); 256 list.push_back(p6);
257 } else if (rdm * totalpnbar < Kinematics 257 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
258 KinematicsUt 258 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
259 KinematicsUt 259 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
260 KinematicsUt 260 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab)) {
261 // Fourth condition 261 // Fourth condition
262 Particle* p1 = new Particle(PiPlus, 262 Particle* p1 = new Particle(PiPlus, zero, rcol);
263 Particle* p2 = new Particle(PiMinus, 263 Particle* p2 = new Particle(PiMinus, zero, rcol);
264 Particle* p3 = new Particle(PiZero, 264 Particle* p3 = new Particle(PiZero, zero, rcol);
265 Particle* p4 = new Particle(PiMinus, 265 Particle* p4 = new Particle(PiMinus, zero, rcol);
266 266
267 list.push_back(p1); 267 list.push_back(p1);
268 list.push_back(p2); 268 list.push_back(p2);
269 list.push_back(p3); 269 list.push_back(p3);
270 list.push_back(p4); 270 list.push_back(p4);
271 } else if (rdm * totalpnbar < Kinematics 271 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
272 KinematicsUt 272 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
273 KinematicsUt 273 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
274 KinematicsUt 274 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
275 KinematicsUt 275 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab)) {
276 // Fifth condition 276 // Fifth condition
277 Particle* p1 = new Particle(PiPlus, 277 Particle* p1 = new Particle(PiPlus, zero, rcol);
278 Particle* p2 = new Particle(PiMinus, 278 Particle* p2 = new Particle(PiMinus, zero, rcol);
279 Particle* p3 = new Particle(PiZero, 279 Particle* p3 = new Particle(PiZero, zero, rcol);
280 Particle* p4 = new Particle(KMinus, 280 Particle* p4 = new Particle(KMinus, zero, rcol);
281 Particle* p5 = new Particle(KZero, z 281 Particle* p5 = new Particle(KZero, zero, rcol);
282 282
283 list.push_back(p1); 283 list.push_back(p1);
284 list.push_back(p2); 284 list.push_back(p2);
285 list.push_back(p3); 285 list.push_back(p3);
286 list.push_back(p4); 286 list.push_back(p4);
287 list.push_back(p5); 287 list.push_back(p5);
288 } else if (rdm * totalpnbar < Kinematics 288 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
289 KinematicsUt 289 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
290 KinematicsUt 290 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
291 KinematicsUt 291 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
292 KinematicsUt 292 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
293 KinematicsUt 293 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab)) {
294 // Sixth condition 294 // Sixth condition
295 Particle* p1 = new Particle(PiPlus, 295 Particle* p1 = new Particle(PiPlus, zero, rcol);
296 Particle* p2 = new Particle(PiMinus, 296 Particle* p2 = new Particle(PiMinus, zero, rcol);
297 Particle* p3 = new Particle(KMinus, 297 Particle* p3 = new Particle(KMinus, zero, rcol);
298 Particle* p4 = new Particle(KZero, z 298 Particle* p4 = new Particle(KZero, zero, rcol);
299 299
300 list.push_back(p1); 300 list.push_back(p1);
301 list.push_back(p2); 301 list.push_back(p2);
302 list.push_back(p3); 302 list.push_back(p3);
303 list.push_back(p4); 303 list.push_back(p4);
304 } else if (rdm * totalpnbar < Kinematics 304 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
305 KinematicsUt 305 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
306 KinematicsUt 306 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
307 KinematicsUt 307 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
308 KinematicsUt 308 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
309 KinematicsUt 309 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
310 KinematicsUt 310 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab)) {
311 // Seventh condition 311 // Seventh condition
312 Particle* p1 = new Particle(PiPlus, 312 Particle* p1 = new Particle(PiPlus, zero, rcol);
313 Particle* p2 = new Particle(PiPlus, 313 Particle* p2 = new Particle(PiPlus, zero, rcol);
314 Particle* p3 = new Particle(PiMinus, 314 Particle* p3 = new Particle(PiMinus, zero, rcol);
315 Particle* p4 = new Particle(PiMinus, 315 Particle* p4 = new Particle(PiMinus, zero, rcol);
316 Particle* p5 = new Particle(PiMinus, 316 Particle* p5 = new Particle(PiMinus, zero, rcol);
317 Particle* p6 = new Particle(PiZero, 317 Particle* p6 = new Particle(PiZero, zero, rcol);
318 318
319 list.push_back(p1); 319 list.push_back(p1);
320 list.push_back(p2); 320 list.push_back(p2);
321 list.push_back(p3); 321 list.push_back(p3);
322 list.push_back(p4); 322 list.push_back(p4);
323 list.push_back(p5); 323 list.push_back(p5);
324 list.push_back(p6); 324 list.push_back(p6);
325 } else if (rdm * totalpnbar < Kinematics 325 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
326 KinematicsUt 326 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
327 KinematicsUt 327 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
328 KinematicsUt 328 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
329 KinematicsUt 329 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
330 KinematicsUt 330 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
331 KinematicsUt 331 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab) +
332 KinematicsUt 332 KinematicsUtils::compute_xs(NPbar_2pip_3pim, plab)) {
333 // Eighth condition 333 // Eighth condition
334 Particle* p1 = new Particle(PiPlus, 334 Particle* p1 = new Particle(PiPlus, zero, rcol);
335 Particle* p2 = new Particle(PiPlus, 335 Particle* p2 = new Particle(PiPlus, zero, rcol);
336 Particle* p3 = new Particle(PiMinus, 336 Particle* p3 = new Particle(PiMinus, zero, rcol);
337 Particle* p4 = new Particle(PiMinus, 337 Particle* p4 = new Particle(PiMinus, zero, rcol);
338 Particle* p5 = new Particle(PiMinus, 338 Particle* p5 = new Particle(PiMinus, zero, rcol);
339 339
340 list.push_back(p1); 340 list.push_back(p1);
341 list.push_back(p2); 341 list.push_back(p2);
342 list.push_back(p3); 342 list.push_back(p3);
343 list.push_back(p4); 343 list.push_back(p4);
344 list.push_back(p5); 344 list.push_back(p5);
345 } else { 345 } else {
346 // Default condition 346 // Default condition
347 //std::cout << "default condition pnba 347 //std::cout << "default condition pnbar"<< std::endl;
348 if(rdm < (1.-kaonicFSprob)){ // pionic 348 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
349 INCL_DEBUG("pionic npbar final s 349 INCL_DEBUG("pionic npbar final state chosen" << '\n');
350 sum = read_file(dataPathnpbar, p 350 sum = read_file(dataPathnpbar, probabilities, particle_types);
351 rdm = (rdm/(1.-kaonicFSprob))*su 351 rdm = (rdm/(1.-kaonicFSprob))*sum; //normalize by the sum of probabilities in the file
352 //now get the line number in the 352 //now get the line number in the file where the FS particles are stored:
353 G4int n = findStringNumber(rdm, 353 G4int n = findStringNumber(rdm, probabilities)-1;
354 for(G4int j = 0; j < static_cast 354 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
355 if(particle_types[n][j] == "pi 355 if(particle_types[n][j] == "pi0"){
356 Particle *p = new Particle(P 356 Particle *p = new Particle(PiZero, zero, rcol);
357 list.push_back(p); 357 list.push_back(p);
358 } 358 }
359 else if(particle_types[n][j] = 359 else if(particle_types[n][j] == "pi-"){
360 Particle *p = new Particle(P 360 Particle *p = new Particle(PiMinus, zero, rcol);
361 list.push_back(p); 361 list.push_back(p);
362 } 362 }
363 else if(particle_types[n][j] = 363 else if(particle_types[n][j] == "pi+"){
364 Particle *p = new Particle(P 364 Particle *p = new Particle(PiPlus, zero, rcol);
365 list.push_back(p); 365 list.push_back(p);
366 } 366 }
367 else if(particle_types[n][j] = 367 else if(particle_types[n][j] == "omega"){
368 Particle *p = new Particle(O 368 Particle *p = new Particle(Omega, zero, rcol);
369 list.push_back(p); 369 list.push_back(p);
370 } 370 }
371 else if(particle_types[n][j] = 371 else if(particle_types[n][j] == "eta"){
372 Particle *p = new Particle(E 372 Particle *p = new Particle(Eta, zero, rcol);
373 list.push_back(p); 373 list.push_back(p);
374 } 374 }
375 else{ 375 else{
376 INCL_ERROR("Some non-existin 376 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
377 for(G4int jj = 0; jj < stati 377 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
378 std::cout << "gotcha! " << 378 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
379 } 379 }
380 } 380 }
381 } 381 }
382 } // end of pionic option 382 } // end of pionic option
383 else{ 383 else{
384 INCL_DEBUG("kaonic npbar final s 384 INCL_DEBUG("kaonic npbar final state chosen" << '\n');
385 sum = read_file(dataPathnpbark, 385 sum = read_file(dataPathnpbark, probabilities, particle_types);
386 rdm = ((1-rdm)/kaonicFSprob)*sum 386 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
387 //now get the line number in the 387 //now get the line number in the file where the FS particles are stored:
388 G4int n = findStringNumber(rdm, 388 G4int n = findStringNumber(rdm, probabilities)-1;
389 for(G4int j = 0; j < static_cast 389 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
390 if(particle_types[n][j] == "pi 390 if(particle_types[n][j] == "pi0"){
391 Particle *p = new Particle 391 Particle *p = new Particle(PiZero, zero, rcol);
392 list.push_back(p); 392 list.push_back(p);
393 } 393 }
394 else if(particle_types[n][j] = 394 else if(particle_types[n][j] == "pi-"){
395 Particle *p = new Particle 395 Particle *p = new Particle(PiMinus, zero, rcol);
396 list.push_back(p); 396 list.push_back(p);
397 } 397 }
398 else if(particle_types[n][j] = 398 else if(particle_types[n][j] == "pi+"){
399 Particle *p = new Particle 399 Particle *p = new Particle(PiPlus, zero, rcol);
400 list.push_back(p); 400 list.push_back(p);
401 } 401 }
402 else if(particle_types[n][j] = 402 else if(particle_types[n][j] == "omega"){
403 Particle *p = new Particle 403 Particle *p = new Particle(Omega, zero, rcol);
404 list.push_back(p); 404 list.push_back(p);
405 } 405 }
406 else if(particle_types[n][j] = 406 else if(particle_types[n][j] == "eta"){
407 Particle *p = new Particle 407 Particle *p = new Particle(Eta, zero, rcol);
408 list.push_back(p); 408 list.push_back(p);
409 } 409 }
410 else if(particle_types[n][j] = 410 else if(particle_types[n][j] == "K-"){
411 Particle *p = new Particle 411 Particle *p = new Particle(KMinus, zero, rcol);
412 list.push_back(p); 412 list.push_back(p);
413 } 413 }
414 else if(particle_types[n][j] = 414 else if(particle_types[n][j] == "K+"){
415 Particle *p = new Particle 415 Particle *p = new Particle(KPlus, zero, rcol);
416 list.push_back(p); 416 list.push_back(p);
417 } 417 }
418 else if(particle_types[n][j] = 418 else if(particle_types[n][j] == "K0"){
419 Particle *p = new Particle 419 Particle *p = new Particle(KZero, zero, rcol);
420 list.push_back(p); 420 list.push_back(p);
421 } 421 }
422 else if(particle_types[n][j] = 422 else if(particle_types[n][j] == "K0b"){
423 Particle *p = new Particle 423 Particle *p = new Particle(KZeroBar, zero, rcol);
424 list.push_back(p); 424 list.push_back(p);
425 } 425 }
426 else{ 426 else{
427 INCL_ERROR("Some non-exist 427 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
428 for(G4int jj = 0; jj < sta 428 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
429 std::cout << "gotcha! " 429 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
430 } 430 }
431 } 431 }
432 } 432 }
433 } // end of kaonic option 433 } // end of kaonic option
434 } // end of default annihilation 434 } // end of default annihilation
435 435
436 } 436 }
437 else if(nucleon->getType()==Proton && anti 437 else if(nucleon->getType()==Proton && antinucleon->getType()==antiNeutron){
438 const G4double totalpnbar = KinematicsUt 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 439 // xs is same for npbar, but the fs has different charge
440 440
441 if (rdm * totalpnbar < KinematicsUtils:: 441 if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab)) {
442 // First condition 442 // First condition
443 Particle* p1 = new Particle(PiPlus, 443 Particle* p1 = new Particle(PiPlus, zero, rcol);
444 Particle* p2 = new Particle(PiMinus, 444 Particle* p2 = new Particle(PiMinus, zero, rcol);
445 Particle* p3 = new Particle(PiPlus, 445 Particle* p3 = new Particle(PiPlus, zero, rcol);
446 446
447 list.push_back(p1); 447 list.push_back(p1);
448 list.push_back(p2); 448 list.push_back(p2);
449 list.push_back(p3); 449 list.push_back(p3);
450 } else if (rdm * totalpnbar < Kinematics 450 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
451 KinematicsUt 451 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab)) {
452 // Second condition 452 // Second condition
453 Particle* p1 = new Particle(PiPlus, 453 Particle* p1 = new Particle(PiPlus, zero, rcol);
454 Particle* p2 = new Particle(PiMinus, 454 Particle* p2 = new Particle(PiMinus, zero, rcol);
455 Particle* p3 = new Particle(PiZero, 455 Particle* p3 = new Particle(PiZero, zero, rcol);
456 Particle* p4 = new Particle(PiZero, 456 Particle* p4 = new Particle(PiZero, zero, rcol);
457 Particle* p5 = new Particle(PiPlus, 457 Particle* p5 = new Particle(PiPlus, zero, rcol);
458 458
459 list.push_back(p1); 459 list.push_back(p1);
460 list.push_back(p2); 460 list.push_back(p2);
461 list.push_back(p3); 461 list.push_back(p3);
462 list.push_back(p4); 462 list.push_back(p4);
463 list.push_back(p5); 463 list.push_back(p5);
464 } else if (rdm * totalpnbar < Kinematics 464 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
465 KinematicsUt 465 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
466 KinematicsUt 466 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab)) {
467 // Third condition 467 // Third condition
468 Particle* p1 = new Particle(PiPlus, 468 Particle* p1 = new Particle(PiPlus, zero, rcol);
469 Particle* p2 = new Particle(PiMinus, 469 Particle* p2 = new Particle(PiMinus, zero, rcol);
470 Particle* p3 = new Particle(PiZero, 470 Particle* p3 = new Particle(PiZero, zero, rcol);
471 Particle* p4 = new Particle(PiZero, 471 Particle* p4 = new Particle(PiZero, zero, rcol);
472 Particle* p5 = new Particle(PiZero, 472 Particle* p5 = new Particle(PiZero, zero, rcol);
473 Particle* p6 = new Particle(PiPlus, 473 Particle* p6 = new Particle(PiPlus, zero, rcol);
474 474
475 list.push_back(p1); 475 list.push_back(p1);
476 list.push_back(p2); 476 list.push_back(p2);
477 list.push_back(p3); 477 list.push_back(p3);
478 list.push_back(p4); 478 list.push_back(p4);
479 list.push_back(p5); 479 list.push_back(p5);
480 list.push_back(p6); 480 list.push_back(p6);
481 } else if (rdm * totalpnbar < Kinematics 481 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
482 KinematicsUt 482 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
483 KinematicsUt 483 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
484 KinematicsUt 484 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab)) {
485 // Fourth condition 485 // Fourth condition
486 Particle* p1 = new Particle(PiPlus, 486 Particle* p1 = new Particle(PiPlus, zero, rcol);
487 Particle* p2 = new Particle(PiMinus, 487 Particle* p2 = new Particle(PiMinus, zero, rcol);
488 Particle* p3 = new Particle(PiZero, 488 Particle* p3 = new Particle(PiZero, zero, rcol);
489 Particle* p4 = new Particle(PiPlus, 489 Particle* p4 = new Particle(PiPlus, zero, rcol);
490 490
491 list.push_back(p1); 491 list.push_back(p1);
492 list.push_back(p2); 492 list.push_back(p2);
493 list.push_back(p3); 493 list.push_back(p3);
494 list.push_back(p4); 494 list.push_back(p4);
495 } else if (rdm * totalpnbar < Kinematics 495 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
496 KinematicsUt 496 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
497 KinematicsUt 497 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
498 KinematicsUt 498 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
499 KinematicsUt 499 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab)) {
500 // Fifth condition 500 // Fifth condition
501 Particle* p1 = new Particle(PiPlus, 501 Particle* p1 = new Particle(PiPlus, zero, rcol);
502 Particle* p2 = new Particle(PiMinus, 502 Particle* p2 = new Particle(PiMinus, zero, rcol);
503 Particle* p3 = new Particle(PiZero, 503 Particle* p3 = new Particle(PiZero, zero, rcol);
504 Particle* p4 = new Particle(KPlus, z 504 Particle* p4 = new Particle(KPlus, zero, rcol);
505 Particle* p5 = new Particle(KZeroBar 505 Particle* p5 = new Particle(KZeroBar, zero, rcol);
506 506
507 list.push_back(p1); 507 list.push_back(p1);
508 list.push_back(p2); 508 list.push_back(p2);
509 list.push_back(p3); 509 list.push_back(p3);
510 list.push_back(p4); 510 list.push_back(p4);
511 list.push_back(p5); 511 list.push_back(p5);
512 } else if (rdm * totalpnbar < Kinematics 512 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
513 KinematicsUt 513 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
514 KinematicsUt 514 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
515 KinematicsUt 515 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
516 KinematicsUt 516 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
517 KinematicsUt 517 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab)) {
518 // Sixth condition 518 // Sixth condition
519 Particle* p1 = new Particle(PiPlus, 519 Particle* p1 = new Particle(PiPlus, zero, rcol);
520 Particle* p2 = new Particle(PiMinus, 520 Particle* p2 = new Particle(PiMinus, zero, rcol);
521 Particle* p3 = new Particle(KPlus, z 521 Particle* p3 = new Particle(KPlus, zero, rcol);
522 Particle* p4 = new Particle(KZeroBar 522 Particle* p4 = new Particle(KZeroBar, zero, rcol);
523 523
524 list.push_back(p1); 524 list.push_back(p1);
525 list.push_back(p2); 525 list.push_back(p2);
526 list.push_back(p3); 526 list.push_back(p3);
527 list.push_back(p4); 527 list.push_back(p4);
528 } else if (rdm * totalpnbar < Kinematics 528 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
529 KinematicsUt 529 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
530 KinematicsUt 530 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
531 KinematicsUt 531 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
532 KinematicsUt 532 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
533 KinematicsUt 533 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
534 KinematicsUt 534 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab)) {
535 // Seventh condition 535 // Seventh condition
536 Particle* p1 = new Particle(PiPlus, 536 Particle* p1 = new Particle(PiPlus, zero, rcol);
537 Particle* p2 = new Particle(PiPlus, 537 Particle* p2 = new Particle(PiPlus, zero, rcol);
538 Particle* p3 = new Particle(PiMinus, 538 Particle* p3 = new Particle(PiMinus, zero, rcol);
539 Particle* p4 = new Particle(PiMinus, 539 Particle* p4 = new Particle(PiMinus, zero, rcol);
540 Particle* p5 = new Particle(PiPlus, 540 Particle* p5 = new Particle(PiPlus, zero, rcol);
541 Particle* p6 = new Particle(PiZero, 541 Particle* p6 = new Particle(PiZero, zero, rcol);
542 542
543 list.push_back(p1); 543 list.push_back(p1);
544 list.push_back(p2); 544 list.push_back(p2);
545 list.push_back(p3); 545 list.push_back(p3);
546 list.push_back(p4); 546 list.push_back(p4);
547 list.push_back(p5); 547 list.push_back(p5);
548 list.push_back(p6); 548 list.push_back(p6);
549 } else if (rdm * totalpnbar < Kinematics 549 } else if (rdm * totalpnbar < KinematicsUtils::compute_xs(NPbar_pip_2pim, plab) +
550 KinematicsUt 550 KinematicsUtils::compute_xs(NPbar_pip_2pim_2pi0, plab) +
551 KinematicsUt 551 KinematicsUtils::compute_xs(NPbar_pip_2pim_3pi0, plab) +
552 KinematicsUt 552 KinematicsUtils::compute_xs(NPbar_pip_2pim_pi0, plab) +
553 KinematicsUt 553 KinematicsUtils::compute_xs(NPbar_pip_pim_pi0_Km_K0, plab) +
554 KinematicsUt 554 KinematicsUtils::compute_xs(NPbar_pip_pim_Km_K0, plab) +
555 KinematicsUt 555 KinematicsUtils::compute_xs(NPbar_2pip_3pim_pi0, plab) +
556 KinematicsUt 556 KinematicsUtils::compute_xs(NPbar_2pip_3pim, plab)) {
557 // Eighth condition 557 // Eighth condition
558 Particle* p1 = new Particle(PiPlus, 558 Particle* p1 = new Particle(PiPlus, zero, rcol);
559 Particle* p2 = new Particle(PiPlus, 559 Particle* p2 = new Particle(PiPlus, zero, rcol);
560 Particle* p3 = new Particle(PiMinus, 560 Particle* p3 = new Particle(PiMinus, zero, rcol);
561 Particle* p4 = new Particle(PiMinus, 561 Particle* p4 = new Particle(PiMinus, zero, rcol);
562 Particle* p5 = new Particle(PiPlus, 562 Particle* p5 = new Particle(PiPlus, zero, rcol);
563 563
564 list.push_back(p1); 564 list.push_back(p1);
565 list.push_back(p2); 565 list.push_back(p2);
566 list.push_back(p3); 566 list.push_back(p3);
567 list.push_back(p4); 567 list.push_back(p4);
568 list.push_back(p5); 568 list.push_back(p5);
569 } else { 569 } else {
570 // Default condition 570 // Default condition
571 if(rdm < (1.-kaonicFSprob)){ // pionic 571 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
572 INCL_DEBUG("pionic pnbar final s 572 INCL_DEBUG("pionic pnbar final state chosen" << '\n');
573 sum = read_file(dataPathpnbar, p 573 sum = read_file(dataPathpnbar, probabilities, particle_types);
574 rdm = (rdm/(1.-kaonicFSprob))*su 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 575 //now get the line number in the file where the FS particles are stored:
576 G4int n = findStringNumber(rdm, 576 G4int n = findStringNumber(rdm, probabilities)-1;
577 for(G4int j = 0; j < static_cast 577 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
578 if(particle_types[n][j] == "pi 578 if(particle_types[n][j] == "pi0"){
579 Particle *p = new Particle(P 579 Particle *p = new Particle(PiZero, zero, rcol);
580 list.push_back(p); 580 list.push_back(p);
581 } 581 }
582 else if(particle_types[n][j] = 582 else if(particle_types[n][j] == "pi-"){
583 Particle *p = new Particle(P 583 Particle *p = new Particle(PiMinus, zero, rcol);
584 list.push_back(p); 584 list.push_back(p);
585 } 585 }
586 else if(particle_types[n][j] = 586 else if(particle_types[n][j] == "pi+"){
587 Particle *p = new Particle(P 587 Particle *p = new Particle(PiPlus, zero, rcol);
588 list.push_back(p); 588 list.push_back(p);
589 } 589 }
590 else if(particle_types[n][j] = 590 else if(particle_types[n][j] == "omega"){
591 Particle *p = new Particle(O 591 Particle *p = new Particle(Omega, zero, rcol);
592 list.push_back(p); 592 list.push_back(p);
593 } 593 }
594 else if(particle_types[n][j] = 594 else if(particle_types[n][j] == "eta"){
595 Particle *p = new Particle(E 595 Particle *p = new Particle(Eta, zero, rcol);
596 list.push_back(p); 596 list.push_back(p);
597 } 597 }
598 else{ 598 else{
599 INCL_ERROR("Some non-existin 599 INCL_ERROR("Some non-existing FS particle detected when reading nbar FS files");
600 for(G4int jj = 0; jj < stati 600 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
601 std::cout << "gotcha! " << 601 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
602 } 602 }
603 } 603 }
604 } 604 }
605 } // end of pionic option 605 } // end of pionic option
606 else{ 606 else{
607 INCL_DEBUG("kaonic pnbar final s 607 INCL_DEBUG("kaonic pnbar final state chosen" << '\n');
608 sum = read_file(dataPathnpbark, 608 sum = read_file(dataPathnpbark, probabilities, particle_types);
609 rdm = ((1-rdm)/kaonicFSprob)*sum 609 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
610 //now get the line number in the 610 //now get the line number in the file where the FS particles are stored:
611 G4int n = findStringNumber(rdm, 611 G4int n = findStringNumber(rdm, probabilities)-1;
612 for(G4int j = 0; j < static_cast 612 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
613 if(particle_types[n][j] == "pi 613 if(particle_types[n][j] == "pi0"){
614 Particle *p = new Particle 614 Particle *p = new Particle(PiZero, zero, rcol);
615 list.push_back(p); 615 list.push_back(p);
616 } 616 }
617 else if(particle_types[n][j] = 617 else if(particle_types[n][j] == "pi-"){
618 Particle *p = new Particle 618 Particle *p = new Particle(PiMinus, zero, rcol);
619 list.push_back(p); 619 list.push_back(p);
620 } 620 }
621 else if(particle_types[n][j] = 621 else if(particle_types[n][j] == "pi+"){
622 Particle *p = new Particle 622 Particle *p = new Particle(PiPlus, zero, rcol);
623 list.push_back(p); 623 list.push_back(p);
624 } 624 }
625 else if(particle_types[n][j] = 625 else if(particle_types[n][j] == "omega"){
626 Particle *p = new Particle 626 Particle *p = new Particle(Omega, zero, rcol);
627 list.push_back(p); 627 list.push_back(p);
628 } 628 }
629 else if(particle_types[n][j] = 629 else if(particle_types[n][j] == "eta"){
630 Particle *p = new Particle 630 Particle *p = new Particle(Eta, zero, rcol);
631 list.push_back(p); 631 list.push_back(p);
632 } 632 }
633 else if(particle_types[n][j] = 633 else if(particle_types[n][j] == "K-"){
634 Particle *p = new Particle 634 Particle *p = new Particle(KMinus, zero, rcol);
635 list.push_back(p); 635 list.push_back(p);
636 } 636 }
637 else if(particle_types[n][j] = 637 else if(particle_types[n][j] == "K+"){
638 Particle *p = new Particle 638 Particle *p = new Particle(KPlus, zero, rcol);
639 list.push_back(p); 639 list.push_back(p);
640 } 640 }
641 else if(particle_types[n][j] = 641 else if(particle_types[n][j] == "K0"){
642 Particle *p = new Particle 642 Particle *p = new Particle(KZero, zero, rcol);
643 list.push_back(p); 643 list.push_back(p);
644 } 644 }
645 else if(particle_types[n][j] = 645 else if(particle_types[n][j] == "K0b"){
646 Particle *p = new Particle 646 Particle *p = new Particle(KZeroBar, zero, rcol);
647 list.push_back(p); 647 list.push_back(p);
648 } 648 }
649 else{ 649 else{
650 INCL_ERROR("Some non-exist 650 INCL_ERROR("Some non-existing FS particle detected when reading pnbar FS files");
651 for(G4int jj = 0; jj < sta 651 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
652 std::cout << "gotcha! " 652 std::cout << "gotcha! " << particle_types[n][jj] << std::endl; }
653 } 653 }
654 } 654 }
655 } // end of kaonic option 655 } // end of kaonic option
656 } // end of default annihilation 656 } // end of default annihilation
657 657
658 } 658 }
659 else{ //ppbar or nnbar 659 else{ //ppbar or nnbar
660 //std::cout << "ppbar or nnbar"<< std::e 660 //std::cout << "ppbar or nnbar"<< std::endl;
661 const G4double totalppbar = KinematicsUt 661 const G4double totalppbar = KinematicsUtils::compute_xs(BFMM6, plab);
662 // same for nnbar 662 // same for nnbar
663 663
664 if (rdm * totalppbar < KinematicsUtils:: 664 if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab)) {
665 // First condition 665 // First condition
666 Particle* p1 = new Particle(PiPlus, 666 Particle* p1 = new Particle(PiPlus, zero, rcol);
667 Particle* p2 = new Particle(PiMinus, 667 Particle* p2 = new Particle(PiMinus, zero, rcol);
668 668
669 list.push_back(p1); 669 list.push_back(p1);
670 list.push_back(p2); 670 list.push_back(p2);
671 671
672 672
673 } else if (rdm * totalppbar < Kinematics 673 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
674 KinematicsUt 674 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab)) {
675 // Second condition 675 // Second condition
676 Particle* p1 = new Particle(PiPlus, 676 Particle* p1 = new Particle(PiPlus, zero, rcol);
677 Particle* p2 = new Particle(PiMinus, 677 Particle* p2 = new Particle(PiMinus, zero, rcol);
678 Particle* p3 = new Particle(PiZero, 678 Particle* p3 = new Particle(PiZero, zero, rcol);
679 679
680 list.push_back(p1); 680 list.push_back(p1);
681 list.push_back(p2); 681 list.push_back(p2);
682 list.push_back(p3); 682 list.push_back(p3);
683 } else if (rdm * totalppbar < Kinematics 683 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
684 KinematicsUt 684 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
685 KinematicsUt 685 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab)) {
686 // Third condition 686 // Third condition
687 Particle* p1 = new Particle(PiPlus, 687 Particle* p1 = new Particle(PiPlus, zero, rcol);
688 Particle* p2 = new Particle(PiMinus, 688 Particle* p2 = new Particle(PiMinus, zero, rcol);
689 Particle* p3 = new Particle(Omega, z 689 Particle* p3 = new Particle(Omega, zero, rcol);
690 690
691 list.push_back(p1); 691 list.push_back(p1);
692 list.push_back(p2); 692 list.push_back(p2);
693 list.push_back(p3); 693 list.push_back(p3);
694 } else if (rdm * totalppbar < Kinematics 694 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
695 KinematicsUt 695 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
696 KinematicsUt 696 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
697 KinematicsUt 697 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab)) {
698 // Fourth condition 698 // Fourth condition
699 Particle* p1 = new Particle(PiPlus, 699 Particle* p1 = new Particle(PiPlus, zero, rcol);
700 Particle* p2 = new Particle(PiMinus, 700 Particle* p2 = new Particle(PiMinus, zero, rcol);
701 Particle* p3 = new Particle(KPlus, z 701 Particle* p3 = new Particle(KPlus, zero, rcol);
702 Particle* p4 = new Particle(KMinus, 702 Particle* p4 = new Particle(KMinus, zero, rcol);
703 703
704 list.push_back(p1); 704 list.push_back(p1);
705 list.push_back(p2); 705 list.push_back(p2);
706 list.push_back(p3); 706 list.push_back(p3);
707 list.push_back(p4); 707 list.push_back(p4);
708 } else if (rdm * totalppbar < Kinematics 708 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
709 KinematicsUt 709 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
710 KinematicsUt 710 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
711 KinematicsUt 711 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
712 KinematicsUt 712 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab)) {
713 // Fifth condition 713 // Fifth condition
714 Particle* p1 = new Particle(PiPlus, 714 Particle* p1 = new Particle(PiPlus, zero, rcol);
715 Particle* p2 = new Particle(PiMinus, 715 Particle* p2 = new Particle(PiMinus, zero, rcol);
716 Particle* p3 = new Particle(PiZero, 716 Particle* p3 = new Particle(PiZero, zero, rcol);
717 Particle* p4 = new Particle(KPlus, z 717 Particle* p4 = new Particle(KPlus, zero, rcol);
718 Particle* p5 = new Particle(KMinus, 718 Particle* p5 = new Particle(KMinus, zero, rcol);
719 719
720 list.push_back(p1); 720 list.push_back(p1);
721 list.push_back(p2); 721 list.push_back(p2);
722 list.push_back(p3); 722 list.push_back(p3);
723 list.push_back(p4); 723 list.push_back(p4);
724 list.push_back(p5); 724 list.push_back(p5);
725 } else if (rdm * totalppbar < Kinematics 725 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
726 KinematicsUt 726 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
727 KinematicsUt 727 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
728 KinematicsUt 728 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
729 KinematicsUt 729 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
730 KinematicsUt 730 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab)) {
731 // Sixth condition 731 // Sixth condition
732 Particle* p1 = new Particle(PiPlus, 732 Particle* p1 = new Particle(PiPlus, zero, rcol);
733 Particle* p2 = new Particle(PiMinus, 733 Particle* p2 = new Particle(PiMinus, zero, rcol);
734 Particle* p3 = new Particle(PiPlus, 734 Particle* p3 = new Particle(PiPlus, zero, rcol);
735 Particle* p4 = new Particle(PiMinus, 735 Particle* p4 = new Particle(PiMinus, zero, rcol);
736 736
737 list.push_back(p1); 737 list.push_back(p1);
738 list.push_back(p2); 738 list.push_back(p2);
739 list.push_back(p3); 739 list.push_back(p3);
740 list.push_back(p4); 740 list.push_back(p4);
741 } else if (rdm * totalppbar < Kinematics 741 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
742 KinematicsUt 742 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
743 KinematicsUt 743 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
744 KinematicsUt 744 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
745 KinematicsUt 745 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
746 KinematicsUt 746 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
747 KinematicsUt 747 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab)) {
748 // Seventh condition 748 // Seventh condition
749 Particle* p1 = new Particle(PiPlus, 749 Particle* p1 = new Particle(PiPlus, zero, rcol);
750 Particle* p2 = new Particle(PiMinus, 750 Particle* p2 = new Particle(PiMinus, zero, rcol);
751 Particle* p3 = new Particle(PiPlus, 751 Particle* p3 = new Particle(PiPlus, zero, rcol);
752 Particle* p4 = new Particle(PiMinus, 752 Particle* p4 = new Particle(PiMinus, zero, rcol);
753 Particle* p5 = new Particle(PiZero, 753 Particle* p5 = new Particle(PiZero, zero, rcol);
754 754
755 list.push_back(p1); 755 list.push_back(p1);
756 list.push_back(p2); 756 list.push_back(p2);
757 list.push_back(p3); 757 list.push_back(p3);
758 list.push_back(p4); 758 list.push_back(p4);
759 list.push_back(p5); 759 list.push_back(p5);
760 } else if (rdm * totalppbar < Kinematics 760 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
761 KinematicsUtils::c 761 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
762 KinematicsUtils::c 762 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
763 KinematicsUtils::c 763 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
764 KinematicsUtils::c 764 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
765 KinematicsUtils::c 765 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
766 KinematicsUtils::c 766 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
767 KinematicsUtils::c 767 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab)) {
768 // Eighth condition 768 // Eighth condition
769 Particle* p1 = new Particle(PiPlus, 769 Particle* p1 = new Particle(PiPlus, zero, rcol);
770 Particle* p2 = new Particle(PiMinus, 770 Particle* p2 = new Particle(PiMinus, zero, rcol);
771 Particle* p3 = new Particle(PiPlus, 771 Particle* p3 = new Particle(PiPlus, zero, rcol);
772 Particle* p4 = new Particle(PiMinus, 772 Particle* p4 = new Particle(PiMinus, zero, rcol);
773 Particle* p5 = new Particle(PiZero, 773 Particle* p5 = new Particle(PiZero, zero, rcol);
774 Particle* p6 = new Particle(PiZero, 774 Particle* p6 = new Particle(PiZero, zero, rcol);
775 Particle* p7 = new Particle(PiZero, 775 Particle* p7 = new Particle(PiZero, zero, rcol);
776 776
777 list.push_back(p1); 777 list.push_back(p1);
778 list.push_back(p2); 778 list.push_back(p2);
779 list.push_back(p3); 779 list.push_back(p3);
780 list.push_back(p4); 780 list.push_back(p4);
781 list.push_back(p5); 781 list.push_back(p5);
782 list.push_back(p6); 782 list.push_back(p6);
783 list.push_back(p7); 783 list.push_back(p7);
784 } else if (rdm * totalppbar < Kinematics 784 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
785 KinematicsUt 785 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
786 KinematicsUt 786 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
787 KinematicsUt 787 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
788 KinematicsUt 788 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
789 KinematicsUt 789 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
790 KinematicsUt 790 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
791 KinematicsUt 791 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
792 KinematicsUt 792 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab)) {
793 // Ninth condition 793 // Ninth condition
794 Particle* p1 = new Particle(PiPlus, 794 Particle* p1 = new Particle(PiPlus, zero, rcol);
795 Particle* p2 = new Particle(PiMinus, 795 Particle* p2 = new Particle(PiMinus, zero, rcol);
796 Particle* p3 = new Particle(PiPlus, 796 Particle* p3 = new Particle(PiPlus, zero, rcol);
797 Particle* p4 = new Particle(PiMinus, 797 Particle* p4 = new Particle(PiMinus, zero, rcol);
798 Particle* p5 = new Particle(PiPlus, 798 Particle* p5 = new Particle(PiPlus, zero, rcol);
799 Particle* p6 = new Particle(PiMinus, 799 Particle* p6 = new Particle(PiMinus, zero, rcol);
800 800
801 list.push_back(p1); 801 list.push_back(p1);
802 list.push_back(p2); 802 list.push_back(p2);
803 list.push_back(p3); 803 list.push_back(p3);
804 list.push_back(p4); 804 list.push_back(p4);
805 list.push_back(p5); 805 list.push_back(p5);
806 list.push_back(p6); 806 list.push_back(p6);
807 } else if (rdm * totalppbar < Kinematics 807 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
808 KinematicsUt 808 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
809 KinematicsUt 809 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
810 KinematicsUt 810 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
811 KinematicsUt 811 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
812 KinematicsUt 812 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
813 KinematicsUt 813 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
814 KinematicsUt 814 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
815 KinematicsUt 815 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
816 KinematicsUt 816 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab)) {
817 // Tenth condition 817 // Tenth condition
818 Particle* p1 = new Particle(PiPlus, 818 Particle* p1 = new Particle(PiPlus, zero, rcol);
819 Particle* p2 = new Particle(PiMinus, 819 Particle* p2 = new Particle(PiMinus, zero, rcol);
820 Particle* p3 = new Particle(PiPlus, 820 Particle* p3 = new Particle(PiPlus, zero, rcol);
821 Particle* p4 = new Particle(PiMinus, 821 Particle* p4 = new Particle(PiMinus, zero, rcol);
822 Particle* p5 = new Particle(PiPlus, 822 Particle* p5 = new Particle(PiPlus, zero, rcol);
823 Particle* p6 = new Particle(PiMinus, 823 Particle* p6 = new Particle(PiMinus, zero, rcol);
824 Particle* p7 = new Particle(PiZero, 824 Particle* p7 = new Particle(PiZero, zero, rcol);
825 825
826 list.push_back(p1); 826 list.push_back(p1);
827 list.push_back(p2); 827 list.push_back(p2);
828 list.push_back(p3); 828 list.push_back(p3);
829 list.push_back(p4); 829 list.push_back(p4);
830 list.push_back(p5); 830 list.push_back(p5);
831 list.push_back(p6); 831 list.push_back(p6);
832 list.push_back(p7); 832 list.push_back(p7);
833 } else if (rdm * totalppbar < Kinematics 833 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
834 KinematicsUt 834 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
835 KinematicsUt 835 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
836 KinematicsUt 836 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
837 KinematicsUt 837 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
838 KinematicsUt 838 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
839 KinematicsUt 839 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
840 KinematicsUt 840 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
841 KinematicsUt 841 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
842 KinematicsUt 842 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
843 KinematicsUt 843 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab)) {
844 // Eleventh condition 844 // Eleventh condition
845 Particle* p1 = new Particle(PiPlus, 845 Particle* p1 = new Particle(PiPlus, zero, rcol);
846 Particle* p2 = new Particle(PiMinus, 846 Particle* p2 = new Particle(PiMinus, zero, rcol);
847 Particle* p3 = new Particle(PiPlus, 847 Particle* p3 = new Particle(PiPlus, zero, rcol);
848 Particle* p4 = new Particle(PiMinus, 848 Particle* p4 = new Particle(PiMinus, zero, rcol);
849 Particle* p5 = new Particle(PiPlus, 849 Particle* p5 = new Particle(PiPlus, zero, rcol);
850 Particle* p6 = new Particle(PiMinus, 850 Particle* p6 = new Particle(PiMinus, zero, rcol);
851 Particle* p7 = new Particle(PiZero, 851 Particle* p7 = new Particle(PiZero, zero, rcol);
852 Particle* p8 = new Particle(PiZero, 852 Particle* p8 = new Particle(PiZero, zero, rcol);
853 853
854 list.push_back(p1); 854 list.push_back(p1);
855 list.push_back(p2); 855 list.push_back(p2);
856 list.push_back(p3); 856 list.push_back(p3);
857 list.push_back(p4); 857 list.push_back(p4);
858 list.push_back(p5); 858 list.push_back(p5);
859 list.push_back(p6); 859 list.push_back(p6);
860 list.push_back(p7); 860 list.push_back(p7);
861 list.push_back(p8); 861 list.push_back(p8);
862 } else if (rdm * totalppbar < Kinematics 862 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
863 KinematicsUt 863 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
864 KinematicsUt 864 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
865 KinematicsUt 865 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
866 KinematicsUt 866 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
867 KinematicsUt 867 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
868 KinematicsUt 868 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
869 KinematicsUt 869 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
870 KinematicsUt 870 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
871 KinematicsUt 871 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
872 KinematicsUt 872 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
873 KinematicsUt 873 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab)) {
874 // Twelfth condition 874 // Twelfth condition
875 Particle* p1 = new Particle(PiPlus, 875 Particle* p1 = new Particle(PiPlus, zero, rcol);
876 Particle* p2 = new Particle(PiMinus, 876 Particle* p2 = new Particle(PiMinus, zero, rcol);
877 Particle* p3 = new Particle(PiPlus, 877 Particle* p3 = new Particle(PiPlus, zero, rcol);
878 Particle* p4 = new Particle(PiMinus, 878 Particle* p4 = new Particle(PiMinus, zero, rcol);
879 Particle* p5 = new Particle(PiPlus, 879 Particle* p5 = new Particle(PiPlus, zero, rcol);
880 Particle* p6 = new Particle(PiMinus, 880 Particle* p6 = new Particle(PiMinus, zero, rcol);
881 Particle* p7 = new Particle(PiZero, 881 Particle* p7 = new Particle(PiZero, zero, rcol);
882 Particle* p8 = new Particle(PiZero, 882 Particle* p8 = new Particle(PiZero, zero, rcol);
883 Particle* p9 = new Particle(PiZero, 883 Particle* p9 = new Particle(PiZero, zero, rcol);
884 884
885 list.push_back(p1); 885 list.push_back(p1);
886 list.push_back(p2); 886 list.push_back(p2);
887 list.push_back(p3); 887 list.push_back(p3);
888 list.push_back(p4); 888 list.push_back(p4);
889 list.push_back(p5); 889 list.push_back(p5);
890 list.push_back(p6); 890 list.push_back(p6);
891 list.push_back(p7); 891 list.push_back(p7);
892 list.push_back(p8); 892 list.push_back(p8);
893 list.push_back(p9); 893 list.push_back(p9);
894 } else if (rdm * totalppbar < Kinematics 894 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
895 KinematicsUt 895 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
896 KinematicsUt 896 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
897 KinematicsUt 897 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
898 KinematicsUt 898 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
899 KinematicsUt 899 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
900 KinematicsUt 900 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
901 KinematicsUt 901 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
902 KinematicsUt 902 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
903 KinematicsUt 903 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
904 KinematicsUt 904 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
905 KinematicsUt 905 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab) +
906 KinematicsUt 906 KinematicsUtils::compute_xs(PPbar_4pip_4pim, plab)) {
907 // Thirteenth condition 907 // Thirteenth condition
908 Particle* p1 = new Particle(PiPlus, 908 Particle* p1 = new Particle(PiPlus, zero, rcol);
909 Particle* p2 = new Particle(PiMinus, 909 Particle* p2 = new Particle(PiMinus, zero, rcol);
910 Particle* p3 = new Particle(PiPlus, 910 Particle* p3 = new Particle(PiPlus, zero, rcol);
911 Particle* p4 = new Particle(PiMinus, 911 Particle* p4 = new Particle(PiMinus, zero, rcol);
912 Particle* p5 = new Particle(PiPlus, 912 Particle* p5 = new Particle(PiPlus, zero, rcol);
913 Particle* p6 = new Particle(PiMinus, 913 Particle* p6 = new Particle(PiMinus, zero, rcol);
914 Particle* p7 = new Particle(PiPlus, 914 Particle* p7 = new Particle(PiPlus, zero, rcol);
915 Particle* p8 = new Particle(PiMinus, 915 Particle* p8 = new Particle(PiMinus, zero, rcol);
916 916
917 list.push_back(p1); 917 list.push_back(p1);
918 list.push_back(p2); 918 list.push_back(p2);
919 list.push_back(p3); 919 list.push_back(p3);
920 list.push_back(p4); 920 list.push_back(p4);
921 list.push_back(p5); 921 list.push_back(p5);
922 list.push_back(p6); 922 list.push_back(p6);
923 list.push_back(p7); 923 list.push_back(p7);
924 list.push_back(p8); 924 list.push_back(p8);
925 } else if (rdm * totalppbar < Kinematics 925 } else if (rdm * totalppbar < KinematicsUtils::compute_xs(PPbar_pip_pim, plab) +
926 KinematicsUt 926 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0, plab) +
927 KinematicsUt 927 KinematicsUtils::compute_xs(PPbar_pip_pim_omega, plab) +
928 KinematicsUt 928 KinematicsUtils::compute_xs(PPbar_pip_pim_Kp_Km, plab) +
929 KinematicsUt 929 KinematicsUtils::compute_xs(PPbar_pip_pim_pi0_Kp_Km, plab) +
930 KinematicsUt 930 KinematicsUtils::compute_xs(PPbar_2pip_2pim, plab) +
931 KinematicsUt 931 KinematicsUtils::compute_xs(PPbar_2pip_2pim_pi0, plab) +
932 KinematicsUt 932 KinematicsUtils::compute_xs(PPbar_2pip_2pim_3pi0, plab) +
933 KinematicsUt 933 KinematicsUtils::compute_xs(PPbar_3pip_3pim, plab) +
934 KinematicsUt 934 KinematicsUtils::compute_xs(PPbar_3pip_3pim_pi0, plab) +
935 KinematicsUt 935 KinematicsUtils::compute_xs(PPbar_3pip_3pim_2pi0, plab) +
936 KinematicsUt 936 KinematicsUtils::compute_xs(PPbar_3pip_3pim_3pi0, plab) +
937 KinematicsUt 937 KinematicsUtils::compute_xs(PPbar_4pip_4pim, plab) +
938 KinematicsUt 938 KinematicsUtils::compute_xs(PPbar_4pip_4pim_pi0, plab)) {
939 // Fourteenth condition 939 // Fourteenth condition
940 Particle* p1 = new Particle(PiPlus, 940 Particle* p1 = new Particle(PiPlus, zero, rcol);
941 Particle* p2 = new Particle(PiMinus, 941 Particle* p2 = new Particle(PiMinus, zero, rcol);
942 Particle* p3 = new Particle(PiPlus, 942 Particle* p3 = new Particle(PiPlus, zero, rcol);
943 Particle* p4 = new Particle(PiMinus, 943 Particle* p4 = new Particle(PiMinus, zero, rcol);
944 Particle* p5 = new Particle(PiPlus, 944 Particle* p5 = new Particle(PiPlus, zero, rcol);
945 Particle* p6 = new Particle(PiMinus, 945 Particle* p6 = new Particle(PiMinus, zero, rcol);
946 Particle* p7 = new Particle(PiPlus, 946 Particle* p7 = new Particle(PiPlus, zero, rcol);
947 Particle* p8 = new Particle(PiMinus, 947 Particle* p8 = new Particle(PiMinus, zero, rcol);
948 Particle* p9 = new Particle(PiZero, 948 Particle* p9 = new Particle(PiZero, zero, rcol);
949 949
950 list.push_back(p1); 950 list.push_back(p1);
951 list.push_back(p2); 951 list.push_back(p2);
952 list.push_back(p3); 952 list.push_back(p3);
953 list.push_back(p4); 953 list.push_back(p4);
954 list.push_back(p5); 954 list.push_back(p5);
955 list.push_back(p6); 955 list.push_back(p6);
956 list.push_back(p7); 956 list.push_back(p7);
957 list.push_back(p8); 957 list.push_back(p8);
958 list.push_back(p9); 958 list.push_back(p9);
959 } else { 959 } else {
960 // Default condition 960 // Default condition
961 if(rdm < (1.-kaonicFSprob)){ // pion 961 if(rdm < (1.-kaonicFSprob)){ // pionic FS was chosen
962 INCL_DEBUG("pionic pp final stat 962 INCL_DEBUG("pionic pp final state chosen" << '\n');
963 sum = read_file(dataPathppbar, p 963 sum = read_file(dataPathppbar, probabilities, particle_types);
964 rdm = (rdm/(1.-kaonicFSprob))*su 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 965 //now get the line number in the file where the FS particles are stored:
966 G4int n = findStringNumber(rdm, 966 G4int n = findStringNumber(rdm, probabilities)-1;
967 for(G4int j = 0; j < static_cast 967 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
968 if(particle_types[n][j] == "pi 968 if(particle_types[n][j] == "pi0"){
969 Particle *p = new Particle 969 Particle *p = new Particle(PiZero, zero, rcol);
970 list.push_back(p); 970 list.push_back(p);
971 } 971 }
972 else if(particle_types[n][j] = 972 else if(particle_types[n][j] == "pi-"){
973 Particle *p = new Particle 973 Particle *p = new Particle(PiMinus, zero, rcol);
974 list.push_back(p); 974 list.push_back(p);
975 } 975 }
976 else if(particle_types[n][j] = 976 else if(particle_types[n][j] == "pi+"){
977 Particle *p = new Particle 977 Particle *p = new Particle(PiPlus, zero, rcol);
978 list.push_back(p); 978 list.push_back(p);
979 } 979 }
980 else if(particle_types[n][j] = 980 else if(particle_types[n][j] == "omega"){
981 Particle *p = new Particle 981 Particle *p = new Particle(Omega, zero, rcol);
982 list.push_back(p); 982 list.push_back(p);
983 } 983 }
984 else if(particle_types[n][j] = 984 else if(particle_types[n][j] == "eta"){
985 Particle *p = new Particle 985 Particle *p = new Particle(Eta, zero, rcol);
986 list.push_back(p); 986 list.push_back(p);
987 } 987 }
988 else{ 988 else{
989 INCL_ERROR("Some non-existing FS 989 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
990 for(G4int jj = 0; jj < static_ca 990 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
991 std::cout << "gotcha! " << par 991 std::cout << "gotcha! " << particle_types[n][jj] << std::endl; }
992 } 992 }
993 } 993 }
994 } //end of pionic option 994 } //end of pionic option
995 else{ 995 else{
996 INCL_DEBUG("kaonic pp final stat 996 INCL_DEBUG("kaonic pp final state chosen" << '\n');
997 sum = read_file(dataPathppbark, 997 sum = read_file(dataPathppbark, probabilities, particle_types);
998 rdm = ((1-rdm)/kaonicFSprob)*sum 998 rdm = ((1-rdm)/kaonicFSprob)*sum;//2670 normalize by the sum of probabilities in the file
999 //now get the line number in the 999 //now get the line number in the file where the FS particles are stored:
1000 G4int n = findStringNumber(rdm, 1000 G4int n = findStringNumber(rdm, probabilities)-1;
1001 for(G4int j = 0; j < static_cast< 1001 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
1002 if(particle_types[n][j] == "p 1002 if(particle_types[n][j] == "pi0"){
1003 Particle *p = new Particl 1003 Particle *p = new Particle(PiZero, zero, rcol);
1004 list.push_back(p); 1004 list.push_back(p);
1005 } 1005 }
1006 else if(particle_types[n][j] 1006 else if(particle_types[n][j] == "pi-"){
1007 Particle *p = new Particl 1007 Particle *p = new Particle(PiMinus, zero, rcol);
1008 list.push_back(p); 1008 list.push_back(p);
1009 } 1009 }
1010 else if(particle_types[n][j] 1010 else if(particle_types[n][j] == "pi+"){
1011 Particle *p = new Particl 1011 Particle *p = new Particle(PiPlus, zero, rcol);
1012 list.push_back(p); 1012 list.push_back(p);
1013 } 1013 }
1014 else if(particle_types[n][j] 1014 else if(particle_types[n][j] == "omega"){
1015 Particle *p = new Particl 1015 Particle *p = new Particle(Omega, zero, rcol);
1016 list.push_back(p); 1016 list.push_back(p);
1017 } 1017 }
1018 else if(particle_types[n][j] 1018 else if(particle_types[n][j] == "eta"){
1019 Particle *p = new Particl 1019 Particle *p = new Particle(Eta, zero, rcol);
1020 list.push_back(p); 1020 list.push_back(p);
1021 } 1021 }
1022 else if(particle_types[n][j] 1022 else if(particle_types[n][j] == "K-"){
1023 Particle *p = new Particl 1023 Particle *p = new Particle(KMinus, zero, rcol);
1024 list.push_back(p); 1024 list.push_back(p);
1025 } 1025 }
1026 else if(particle_types[n][j] 1026 else if(particle_types[n][j] == "K+"){
1027 Particle *p = new Particl 1027 Particle *p = new Particle(KPlus, zero, rcol);
1028 list.push_back(p); 1028 list.push_back(p);
1029 } 1029 }
1030 else if(particle_types[n][j] 1030 else if(particle_types[n][j] == "K0"){
1031 Particle *p = new Particl 1031 Particle *p = new Particle(KZero, zero, rcol);
1032 list.push_back(p); 1032 list.push_back(p);
1033 } 1033 }
1034 else if(particle_types[n][j] 1034 else if(particle_types[n][j] == "K0b"){
1035 Particle *p = new Particl 1035 Particle *p = new Particle(KZeroBar, zero, rcol);
1036 list.push_back(p); 1036 list.push_back(p);
1037 } 1037 }
1038 else{ 1038 else{
1039 INCL_ERROR("Some non-exis 1039 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
1040 for(G4int jj = 0; jj < st 1040 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
1041 std::cout << "gotcha! " 1041 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
1042 } 1042 }
1043 } 1043 }
1044 } 1044 }
1045 } // end of kaonic option 1045 } // end of kaonic option
1046 } // end of default condition 1046 } // end of default condition
1047 } // end of ppbar and nnbar case 1047 } // end of ppbar and nnbar case
1048 1048
1049 1049
1050 1050
1051 nucleon->setType(list[0]->getType()); 1051 nucleon->setType(list[0]->getType());
1052 antinucleon->setType(list[1]->getType()); 1052 antinucleon->setType(list[1]->getType());
1053 1053
1054 ParticleList finallist; 1054 ParticleList finallist;
1055 1055
1056 finallist.push_back(nucleon); 1056 finallist.push_back(nucleon);
1057 finallist.push_back(antinucleon); 1057 finallist.push_back(antinucleon);
1058 1058
1059 if(list.size() > 2){ 1059 if(list.size() > 2){
1060 for (G4int i = 2; i < (G4int)(list.size 1060 for (G4int i = 2; i < (G4int)(list.size()); i++) {
1061 finallist.push_back(list[i]); 1061 finallist.push_back(list[i]);
1062 } 1062 }
1063 } 1063 }
1064 1064
1065 if(finallist.size()==2){ 1065 if(finallist.size()==2){
1066 G4double mn=nucleon->getMass(); 1066 G4double mn=nucleon->getMass();
1067 G4double my=antinucleon->getMass(); 1067 G4double my=antinucleon->getMass();
1068 1068
1069 G4double ey=(sqrtS*sqrtS+my*my-mn*mn)/( 1069 G4double ey=(sqrtS*sqrtS+my*my-mn*mn)/(2*sqrtS);
1070 G4double en=std::sqrt(ey*ey-my*my+mn*mn 1070 G4double en=std::sqrt(ey*ey-my*my+mn*mn);
1071 nucleon->setEnergy(en); 1071 nucleon->setEnergy(en);
1072 antinucleon->setEnergy(ey); 1072 antinucleon->setEnergy(ey);
1073 G4double py=std::sqrt(ey*ey-my*my); 1073 G4double py=std::sqrt(ey*ey-my*my);
1074 1074
1075 ThreeVector mom_antinucleon = Random::n 1075 ThreeVector mom_antinucleon = Random::normVector(py);
1076 antinucleon->setMomentum(mom_antinucleo 1076 antinucleon->setMomentum(mom_antinucleon);
1077 nucleon->setMomentum(-mom_antinucleon); 1077 nucleon->setMomentum(-mom_antinucleon);
1078 } 1078 }
1079 else if(finallist.size() > 2){ 1079 else if(finallist.size() > 2){
1080 PhaseSpaceGenerator::generate(sqrtS, fi 1080 PhaseSpaceGenerator::generate(sqrtS, finallist);
1081 } 1081 }
1082 else{ 1082 else{
1083 INCL_ERROR("less than 2 mesons in NNbar 1083 INCL_ERROR("less than 2 mesons in NNbar annihilation!" << '\n');
1084 } 1084 }
1085 1085
1086 1086
1087 for (G4int i = 0; i < 2; i++) { 1087 for (G4int i = 0; i < 2; i++) {
1088 fs->addModifiedParticle(finallist[i]); 1088 fs->addModifiedParticle(finallist[i]);
1089 } 1089 }
1090 if(finallist.size()>2){ 1090 if(finallist.size()>2){
1091 for (G4int i = 2; i < (G4int)(list.size 1091 for (G4int i = 2; i < (G4int)(list.size()); i++) {
1092 fs->addCreatedParticle(finallist[i]); 1092 fs->addCreatedParticle(finallist[i]);
1093 } 1093 }
1094 } 1094 }
1095 1095
1096 1096
1097 //fs->addDestroyedParticle(nucleon); 1097 //fs->addDestroyedParticle(nucleon);
1098 //fs->addDestroyedParticle(antinucleon); 1098 //fs->addDestroyedParticle(antinucleon);
1099 1099
1100 1100
1101 } 1101 }
1102 } 1102 }
1103 1103
1104 1104