| 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 38
39 #include "G4INCLPbarAtrestEntryChannel.hh" 39 #include "G4INCLPbarAtrestEntryChannel.hh"
40 #include "G4INCLRootFinder.hh" 40 #include "G4INCLRootFinder.hh"
41 #include "G4INCLIntersection.hh" 41 #include "G4INCLIntersection.hh"
42 #include "G4INCLCascade.hh" 42 #include "G4INCLCascade.hh"
43 #include <algorithm> 43 #include <algorithm>
44 #include "G4INCLParticle.hh" 44 #include "G4INCLParticle.hh"
45 #include "G4INCLKinematicsUtils.hh" 45 #include "G4INCLKinematicsUtils.hh"
46 #include "G4INCLBinaryCollisionAvatar.hh" 46 #include "G4INCLBinaryCollisionAvatar.hh"
47 #include "G4INCLRandom.hh" 47 #include "G4INCLRandom.hh"
48 #include "G4INCLGlobals.hh" 48 #include "G4INCLGlobals.hh"
49 #include "G4INCLLogger.hh" 49 #include "G4INCLLogger.hh"
50 #include <algorithm> 50 #include <algorithm>
51 #include "G4INCLPhaseSpaceGenerator.hh" 51 #include "G4INCLPhaseSpaceGenerator.hh"
52 #include <iostream> 52 #include <iostream>
53 #include <string> 53 #include <string>
54 #include <sstream> 54 #include <sstream>
55 #include <vector> 55 #include <vector>
56 #include "G4INCLHFB.hh" 56 #include "G4INCLHFB.hh"
57 #include "G4INCLParticleEntryAvatar.hh" 57 #include "G4INCLParticleEntryAvatar.hh"
58 #include "G4INCLNuclearDensityFactory.hh" 58 #include "G4INCLNuclearDensityFactory.hh"
59 #include "G4INCLNDFWoodsSaxon.hh" 59 #include "G4INCLNDFWoodsSaxon.hh"
60 #include "G4INCLNDFModifiedHarmonicOscillator. 60 #include "G4INCLNDFModifiedHarmonicOscillator.hh"
61 #include "G4INCLNDFGaussian.hh" 61 #include "G4INCLNDFGaussian.hh"
62 #include "G4INCLNDFParis.hh" 62 #include "G4INCLNDFParis.hh"
63 #include <string> 63 #include <string>
64 #include <vector> 64 #include <vector>
65 #include <iostream> 65 #include <iostream>
66 #include <fstream> 66 #include <fstream>
67 #include <sstream> 67 #include <sstream>
68 68
69 namespace G4INCL { 69 namespace G4INCL {
70 70
71 PbarAtrestEntryChannel::PbarAtrestEntryChann 71 PbarAtrestEntryChannel::PbarAtrestEntryChannel(Nucleus *n, Particle *p)
72 :theNucleus(n), theParticle(p) 72 :theNucleus(n), theParticle(p)
73 {} 73 {}
74 74
75 PbarAtrestEntryChannel::~PbarAtrestEntryChan 75 PbarAtrestEntryChannel::~PbarAtrestEntryChannel()
76 {} 76 {}
77 77
78 //fill probabilities and particle types from 78 //fill probabilities and particle types from datafile and return probability sum for normalization
79 G4double PbarAtrestEntryChannel::read_file(s 79 G4double PbarAtrestEntryChannel::read_file(std::string filename, std::vector<G4double>& probabilities, std::vector<std::vector<std::string>>& particle_types) {
80 std::ifstream file(filename); 80 std::ifstream file(filename);
81 G4double sum_probs = 0.0; 81 G4double sum_probs = 0.0;
82 if (file.is_open()) { 82 if (file.is_open()) {
83 std::string line; 83 std::string line;
84 while (getline(file, line)) { 84 while (getline(file, line)) {
85 std::istringstream iss(line); 85 std::istringstream iss(line);
86 G4double prob; 86 G4double prob;
87 iss >> prob; 87 iss >> prob;
88 sum_probs += prob; 88 sum_probs += prob;
89 probabilities.push_back(prob); 89 probabilities.push_back(prob);
90 std::vector<std::string> types; 90 std::vector<std::string> types;
91 std::string type; 91 std::string type;
92 while (iss >> type) { 92 while (iss >> type) {
93 types.push_back(type); 93 types.push_back(type);
94 } 94 }
95 particle_types.push_back(std::mo << 95 particle_types.push_back(types);
96 } 96 }
97 } 97 }
98 else std::cout << "ERROR no fread_file " 98 else std::cout << "ERROR no fread_file " << filename << std::endl;
99 99
100 return sum_probs; 100 return sum_probs;
101 } 101 }
102 102
103 //this function will tell you the FS line nu 103 //this function will tell you the FS line number from the datafile based on your random value
104 G4int PbarAtrestEntryChannel::findStringNumb 104 G4int PbarAtrestEntryChannel::findStringNumber(G4double rdm, std::vector<G4double> yields) {
105 G4int stringNumber = -1; 105 G4int stringNumber = -1;
106 G4double smallestsum = 0.0; 106 G4double smallestsum = 0.0;
107 G4double biggestsum = yields[0]; 107 G4double biggestsum = yields[0];
108 //std::cout << "initial input " << rdm << 108 //std::cout << "initial input " << rdm << std::endl;
109 for (G4int i = 0; i < static_cast<G4int>(y 109 for (G4int i = 0; i < static_cast<G4int>(yields.size()-1); i++) {
110 if (rdm >= smallestsum && rdm <= bigge 110 if (rdm >= smallestsum && rdm <= biggestsum) {
111 //std::cout << smallestsum << " an 111 //std::cout << smallestsum << " and " << biggestsum << std::endl;
112 stringNumber = i+1; 112 stringNumber = i+1;
113 } 113 }
114 smallestsum += yields[i]; 114 smallestsum += yields[i];
115 biggestsum += yields[i+1]; 115 biggestsum += yields[i+1];
116 } 116 }
117 if(stringNumber==-1) stringNumber = static 117 if(stringNumber==-1) stringNumber = static_cast<G4int>(yields.size());
118 if(stringNumber==-1){ 118 if(stringNumber==-1){
119 INCL_ERROR("ERROR in findStringNumber (s 119 INCL_ERROR("ERROR in findStringNumber (stringNumber=-1)");
120 std::cout << "ERROR in findStringNumber" 120 std::cout << "ERROR in findStringNumber" << std::endl;
121 } 121 }
122 return stringNumber; 122 return stringNumber;
123 } 123 }
124 124
125 G4double PbarAtrestEntryChannel::fctrl(const 125 G4double PbarAtrestEntryChannel::fctrl(const G4double arg1){
126 G4double factorial=1.0; 126 G4double factorial=1.0;
127 for(G4int k=1; k<=arg1; k++){ 127 for(G4int k=1; k<=arg1; k++){
128 factorial *= k; 128 factorial *= k;
129 } 129 }
130 return factorial; 130 return factorial;
131 } 131 }
132 G4double PbarAtrestEntryChannel::r1(const G4 132 G4double PbarAtrestEntryChannel::r1(const G4int n) {
133 return std::pow(fctrl(2.0*n),-0.5); 133 return std::pow(fctrl(2.0*n),-0.5);
134 } 134 }
135 G4double PbarAtrestEntryChannel::r2(const G4 135 G4double PbarAtrestEntryChannel::r2(const G4int n) {
136 G4int Z = theNucleus->getZ(); 136 G4int Z = theNucleus->getZ();
137 return std::pow((Z/(14.4*n)), 1.5); 137 return std::pow((Z/(14.4*n)), 1.5);
138 } 138 }
139 G4double PbarAtrestEntryChannel::r3(G4double 139 G4double PbarAtrestEntryChannel::r3(G4double x, const G4int n) {
140 G4int Z = theNucleus->getZ(); 140 G4int Z = theNucleus->getZ();
141 return std::pow((x)*Z/(n*14.4),(n-1)); //w 141 return std::pow((x)*Z/(n*14.4),(n-1)); //why x is vector here?
142 } 142 }
143 G4double PbarAtrestEntryChannel::r4(G4double 143 G4double PbarAtrestEntryChannel::r4(G4double x, const G4int n) {
144 G4int Z = theNucleus->getZ(); 144 G4int Z = theNucleus->getZ();
145 return std::exp((-x*Z)/(n*28.8)); 145 return std::exp((-x*Z)/(n*28.8));
146 } 146 }
147 147
148 148
149 G4double PbarAtrestEntryChannel::radial_wave 149 G4double PbarAtrestEntryChannel::radial_wavefunction(G4double x, const G4int n){
150 return r1(n)*r2(n)*r3(x,n)*r4(x,n); //Rad 150 return r1(n)*r2(n)*r3(x,n)*r4(x,n); //Radial wave function par=(n, Z)
151 } 151 }
152 152
153 /* 153 /*
154 G4double PbarAtrestEntryChannel::densityP(G4 154 G4double PbarAtrestEntryChannel::densityP(G4double *x, G4double *par){
155 return 0.16800136/(1.0+std::exp((x[0]-par 155 return 0.16800136/(1.0+std::exp((x[0]-par[2])/par[3])); //P nuclear density
156 */ 156 */
157 G4double PbarAtrestEntryChannel::densityP(G4 157 G4double PbarAtrestEntryChannel::densityP(G4double x) {
158 158
159 const G4bool isProton = ProtonIsTheVic 159 const G4bool isProton = ProtonIsTheVictim();
160 G4int Z = theNucleus->getZ(); //was mo 160 G4int Z = theNucleus->getZ(); //was modified in Cascade.cc
161 G4int A = theNucleus->getA(); //was mo 161 G4int A = theNucleus->getA(); //was modified in Cascade.cc
162 A++; //restoration of original A value 162 A++; //restoration of original A value before annihilation
163 if(isProton == true){Z++;} //restorati 163 if(isProton == true){Z++;} //restoration of original Z value before annihilation
164 164
165 if(A > 19) { 165 if(A > 19) {
166 G4double radius = ParticleTable::get 166 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
167 G4double diffuseness = ParticleTable 167 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Proton, A, Z);
168 G4double maximumRadius = ParticleTab 168 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
169 NuclearDensityFunctions::WoodsSaxon 169 NuclearDensityFunctions::WoodsSaxon rDensityFunction(radius, maximumRadius, diffuseness);
170 return ((x!=0.) ? rDensityFunction(x 170 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
171 } else if(A <= 19 && A > 6) { 171 } else if(A <= 19 && A > 6) {
172 G4double radius = ParticleTable::get 172 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
173 G4double diffuseness = ParticleTable 173 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Proton, A, Z);
174 G4double maximumRadius = ParticleTab 174 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
175 NuclearDensityFunctions::ModifiedHar 175 NuclearDensityFunctions::ModifiedHarmonicOscillator rDensityFunction(radius, maximumRadius, diffuseness);
176 return ((x!=0.) ? rDensityFunction(x 176 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
177 } else if(A <= 6 && A > 2) { // Gaussi 177 } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
178 G4double radius = ParticleTable::get 178 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
179 G4double maximumRadius = ParticleTab 179 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
180 NuclearDensityFunctions::Gaussian rD 180 NuclearDensityFunctions::Gaussian rDensityFunction(maximumRadius, Math::oneOverSqrtThree * radius);
181 return ((x!=0.) ? rDensityFunction(x 181 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
182 } else if(A == 2 && Z == 1) { // densi 182 } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
183 NuclearDensityFunctions::ParisR rDen 183 NuclearDensityFunctions::ParisR rDensityFunction;
184 return ((x!=0.) ? rDensityFunction(x 184 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
185 } else { 185 } else {
186 INCL_ERROR("No nuclear density funct 186 INCL_ERROR("No nuclear density function for target A = "
187 << A << " Z = " << Z << '\n'); 187 << A << " Z = " << Z << '\n');
188 return 0.0; 188 return 0.0;
189 } 189 }
190 190
191 } 191 }
192 /* 192 /*
193 } 193 }
194 G4double PbarAtrestEntryChannel::densityN(G4 194 G4double PbarAtrestEntryChannel::densityN(G4double *x, G4double *par){
195 return 0.16800136/(1.0+std::exp((x[0]-par 195 return 0.16800136/(1.0+std::exp((x[0]-par[2])/par[3])); //N nuclear density
196 } 196 }
197 */ 197 */
198 G4double PbarAtrestEntryChannel::densityN(G4 198 G4double PbarAtrestEntryChannel::densityN(G4double x) {
199 199
200 const G4bool isProton = ProtonIsTheVic 200 const G4bool isProton = ProtonIsTheVictim();
201 G4int Z = theNucleus->getZ(); //was mo 201 G4int Z = theNucleus->getZ(); //was modified in Cascade.cc
202 G4int A = theNucleus->getA(); //was mo 202 G4int A = theNucleus->getA(); //was modified in Cascade.cc
203 A++; //restoration of original A value 203 A++; //restoration of original A value before annihilation
204 if(isProton == true){Z++;} //restorati 204 if(isProton == true){Z++;} //restoration of original Z value before annihilation
205 205
206 if(A > 19) { 206 if(A > 19) {
207 G4double radius = ParticleTable::get 207 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
208 G4double diffuseness = ParticleTable 208 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Neutron, A, Z);
209 G4double maximumRadius = ParticleTab 209 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
210 NuclearDensityFunctions::WoodsSaxon 210 NuclearDensityFunctions::WoodsSaxon rDensityFunction(radius, maximumRadius, diffuseness);
211 return ((x!=0.) ? rDensityFunction(x 211 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
212 } else if(A <= 19 && A > 6) { 212 } else if(A <= 19 && A > 6) {
213 G4double radius = ParticleTable::get 213 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
214 G4double diffuseness = ParticleTable 214 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Neutron, A, Z);
215 G4double maximumRadius = ParticleTab 215 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
216 NuclearDensityFunctions::ModifiedHar 216 NuclearDensityFunctions::ModifiedHarmonicOscillator rDensityFunction(radius, maximumRadius, diffuseness);
217 return ((x!=0.) ? rDensityFunction(x 217 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
218 } else if(A <= 6 && A > 2) { // Gaussi 218 } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
219 G4double radius = ParticleTable::get 219 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
220 G4double maximumRadius = ParticleTab 220 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
221 NuclearDensityFunctions::Gaussian rD 221 NuclearDensityFunctions::Gaussian rDensityFunction(maximumRadius, Math::oneOverSqrtThree * radius);
222 return ((x!=0.) ? rDensityFunction(x 222 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
223 } else if(A == 2 && Z == 1) { // densi 223 } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
224 NuclearDensityFunctions::ParisR rDen 224 NuclearDensityFunctions::ParisR rDensityFunction;
225 return ((x!=0.) ? rDensityFunction(x 225 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
226 } else { 226 } else {
227 INCL_ERROR("No nuclear density funct 227 INCL_ERROR("No nuclear density function for target A = "
228 << A << " Z = " << Z << '\n'); 228 << A << " Z = " << Z << '\n');
229 return 0.0; 229 return 0.0;
230 } 230 }
231 231
232 } 232 }
233 233
234 234
235 G4double PbarAtrestEntryChannel::overlapP(G4 235 G4double PbarAtrestEntryChannel::overlapP(G4double &x, const G4int n){
236 return x*x*r1(n)*r2(n)*r3(x,n)*r4(x,n)*r1 236 return x*x*r1(n)*r2(n)*r3(x,n)*r4(x,n)*r1(n)*r2(n)*r3(x,n)*r4(x,n)*densityP(x);
237 } 237 }
238 G4double PbarAtrestEntryChannel::overlapN(G4 238 G4double PbarAtrestEntryChannel::overlapN(G4double &x, const G4int n){
239 return x*x*r1(n)*r2(n)*r3(x,n)*r4(x,n)*r1 239 return x*x*r1(n)*r2(n)*r3(x,n)*r4(x,n)*r1(n)*r2(n)*r3(x,n)*r4(x,n)*densityN(x);
240 } 240 }
241 241
242 242
243 ParticleList PbarAtrestEntryChannel::makeMes 243 ParticleList PbarAtrestEntryChannel::makeMesonStar() {//This function creates a set of mesons with momenta
244 244
245 // File names 245 // File names
246 #ifdef INCLXX_IN_GEANT4_MODE 246 #ifdef INCLXX_IN_GEANT4_MODE
247 if(!G4FindDataDir("G4INCLDATA")) { 247 if(!G4FindDataDir("G4INCLDATA")) {
248 G4ExceptionDescription ed; 248 G4ExceptionDescription ed;
249 ed << " Data missing: set environment 249 ed << " Data missing: set environment variable G4INCLDATA\n"
250 << " to point to the directory cont 250 << " to point to the directory containing data files needed\n"
251 << " by the INCL++ model" << G4endl 251 << " by the INCL++ model" << G4endl;
252 G4Exception("G4INCLDataFile::readDa 252 G4Exception("G4INCLDataFile::readData()","rawppbarFS.dat, ...",
253 FatalException, ed); 253 FatalException, ed);
254 } 254 }
255 const G4String& dataPath0{G4FindDataDir( << 255 G4String dataPath0{G4FindDataDir("G4INCLDATA")};
256 const G4String& dataPathppbar(dataPath0 << 256 G4String dataPathppbar(dataPath0 + "/rawppbarFS.dat");
257 const G4String& dataPathnpbar(dataPath0 << 257 G4String dataPathnpbar(dataPath0 + "/rawnpbarFS.dat");
258 const G4String& dataPathppbark(dataPath0 << 258 G4String dataPathppbark(dataPath0 + "/rawppbarFSkaonic.dat");
259 const G4String& dataPathnpbark(dataPath0 << 259 G4String dataPathnpbark(dataPath0 + "/rawnpbarFSkaonic.dat");
260 #else 260 #else
261 Config const *theConfig=theNucleus->getS 261 Config const *theConfig=theNucleus->getStore()->getConfig();
262 std::string path; 262 std::string path;
263 if(theConfig) 263 if(theConfig)
264 path = theConfig->getINCLXXDataFilePat 264 path = theConfig->getINCLXXDataFilePath();
265 const std::string& dataPathppbar(path + << 265 std::string dataPathppbar(path + "/rawppbarFS.dat");
266 INCL_DEBUG("Reading https://doi.org/10.1 266 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
267 const std::string& dataPathnpbar(path + << 267 std::string dataPathnpbar(path + "/rawnpbarFS.dat");
268 INCL_DEBUG("Reading https://doi.org/10.1 268 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
269 const std::string& dataPathppbark(path + << 269 std::string dataPathppbark(path + "/rawppbarFSkaonic.dat");
270 INCL_DEBUG("Reading https://doi.org/10.1 270 INCL_DEBUG("Reading https://doi.org/10.1016/j.physrep.2005.03.002 ppbar kaonic final states" << dataPathppbark << '\n');
271 const std::string& dataPathnpbark(path + << 271 std::string dataPathnpbark(path + "/rawnpbarFSkaonic.dat");
272 INCL_DEBUG("Reading https://doi.org/10.1 272 INCL_DEBUG("Reading https://doi.org/10.1007/BF02818764 and https://link.springer.com/article/10.1007/BF02754930 npbar kaonic final states" << dataPathnpbark << '\n');
273 #endif 273 #endif
274 /*std::string path = {"/home/zdemid/INCL/i 274 /*std::string path = {"/home/zdemid/INCL/inclcode/data"};
275 std::string dataPathppbar(path + "/rawppba 275 std::string dataPathppbar(path + "/rawppbarFS.dat");
276 INCL_DEBUG("Reading https://doi.org/10.101 276 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
277 std::string dataPathnpbar(path + "/rawnpba 277 std::string dataPathnpbar(path + "/rawnpbarFS.dat");
278 INCL_DEBUG("Reading https://doi.org/10.101 278 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
279 std::string dataPathppbark(path + "/rawppb 279 std::string dataPathppbark(path + "/rawppbarFSkaonic.dat");
280 INCL_DEBUG("Reading https://doi.org/10.101 280 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar kaonic final states" << dataPathppbark << '\n');
281 std::string dataPathnpbark(path + "/rawnpb 281 std::string dataPathnpbark(path + "/rawnpbarFSkaonic.dat");
282 INCL_DEBUG("Reading https://doi.org/10.101 282 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar kaonic final states" << dataPathnpbark << '\n');
283 */ 283 */
284 //read probabilities and particle types fr 284 //read probabilities and particle types from file
285 std::vector<G4double> probabilities; //wil 285 std::vector<G4double> probabilities; //will store each FS yield
286 std::vector<std::vector<std::string>> part 286 std::vector<std::vector<std::string>> particle_types; //will store particle names
287 G4double sum; //will contain a sum of prob 287 G4double sum; //will contain a sum of probabilities of all FS in the file
288 G4double kaonicFSprob=0.05; //probability 288 G4double kaonicFSprob=0.05; //probability to kave kaonic FS
289 289
290 const G4bool isProton = ProtonIsTheVictim( 290 const G4bool isProton = ProtonIsTheVictim();
291 G4int z = theNucleus->getZ(); //was modifi 291 G4int z = theNucleus->getZ(); //was modified in Cascade.cc
292 G4int a = theNucleus->getA(); //was modifi 292 G4int a = theNucleus->getA(); //was modified in Cascade.cc
293 a++; //restoration of original A value bef 293 a++; //restoration of original A value before annihilation
294 if(isProton == true){z++;} //restoration o 294 if(isProton == true){z++;} //restoration of original Z value before annihilation
295 ThreeVector annihilationPosition; 295 ThreeVector annihilationPosition;
296 ParticleList starlist; 296 ParticleList starlist;
297 ThreeVector mommy; //momentum to be assign 297 ThreeVector mommy; //momentum to be assigned later
298 298
299 //LETS GOOOOOOO!!! 299 //LETS GOOOOOOO!!!
300 G4double rdm = Random::shoot(); 300 G4double rdm = Random::shoot();
301 if(isProton == true){ //protonic annihilat 301 if(isProton == true){ //protonic annihilation
302 INCL_DEBUG("Proton is the victim" << '\n 302 INCL_DEBUG("Proton is the victim" << '\n');
303 if(rdm < (1.-kaonicFSprob)){ // pionic F 303 if(rdm < (1.-kaonicFSprob)){ // pionic FS was chosen
304 INCL_DEBUG("pionic pp final state chos 304 INCL_DEBUG("pionic pp final state chosen" << '\n');
305 sum = read_file(dataPathppbar, probabi 305 sum = read_file(dataPathppbar, probabilities, particle_types);
306 rdm = (rdm/(1.-kaonicFSprob))*sum; //9 306 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.88 normalize by the sum of probabilities in the file
307 //now get the line number in the file 307 //now get the line number in the file where the FS particles are stored:
308 G4int n = findStringNumber(rdm, std::m << 308 G4int n = findStringNumber(rdm, probabilities)-1;
309 if ( n < 0 ) return starlist; 309 if ( n < 0 ) return starlist;
310 for(G4int j = 0; j < static_cast<G4int 310 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
311 if(particle_types[n][j] == "pi0"){ 311 if(particle_types[n][j] == "pi0"){
312 Particle *p = new Particle(PiZero, 312 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
313 starlist.push_back(p); 313 starlist.push_back(p);
314 } 314 }
315 else if(particle_types[n][j] == "pi- 315 else if(particle_types[n][j] == "pi-"){
316 Particle *p = new Particle(PiMinus 316 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
317 starlist.push_back(p); 317 starlist.push_back(p);
318 } 318 }
319 else if(particle_types[n][j] == "pi+ 319 else if(particle_types[n][j] == "pi+"){
320 Particle *p = new Particle(PiPlus, 320 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
321 starlist.push_back(p); 321 starlist.push_back(p);
322 } 322 }
323 else if(particle_types[n][j] == "ome 323 else if(particle_types[n][j] == "omega"){
324 Particle *p = new Particle(Omega, 324 Particle *p = new Particle(Omega, mommy, annihilationPosition);
325 starlist.push_back(p); 325 starlist.push_back(p);
326 } 326 }
327 else if(particle_types[n][j] == "eta 327 else if(particle_types[n][j] == "eta"){
328 Particle *p = new Particle(Eta, mo 328 Particle *p = new Particle(Eta, mommy, annihilationPosition);
329 starlist.push_back(p); 329 starlist.push_back(p);
330 } 330 }
331 else if(particle_types[n][j] == "rho 331 else if(particle_types[n][j] == "rho-"){
332 Particle *p = new Particle(PiMinus 332 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
333 starlist.push_back(p); 333 starlist.push_back(p);
334 Particle *pp = new Particle(PiZero 334 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
335 starlist.push_back(pp); 335 starlist.push_back(pp);
336 } 336 }
337 else if(particle_types[n][j] == "rho 337 else if(particle_types[n][j] == "rho+"){
338 Particle *p = new Particle(PiPlus, 338 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
339 starlist.push_back(p); 339 starlist.push_back(p);
340 Particle *pp = new Particle(PiZero 340 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
341 starlist.push_back(pp); 341 starlist.push_back(pp);
342 } 342 }
343 else if(particle_types[n][j] == "rho 343 else if(particle_types[n][j] == "rho0"){
344 Particle *p = new Particle(PiMinus 344 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
345 starlist.push_back(p); 345 starlist.push_back(p);
346 Particle *pp = new Particle(PiPlus 346 Particle *pp = new Particle(PiPlus, mommy, annihilationPosition);
347 starlist.push_back(pp); 347 starlist.push_back(pp);
348 } 348 }
349 else{ 349 else{
350 INCL_ERROR("Some non-existing FS p 350 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
351 for(G4int jj = 0; jj < static_cast 351 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
352 std::cout << "gotcha! " << parti 352 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
353 } 353 }
354 std::cout << "Some non-existing FS 354 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
355 } 355 }
356 } 356 }
357 } 357 }
358 else{ 358 else{
359 INCL_DEBUG("kaonic pp final state chos 359 INCL_DEBUG("kaonic pp final state chosen" << '\n');
360 sum = read_file(dataPathppbark, probab 360 sum = read_file(dataPathppbark, probabilities, particle_types);
361 rdm = ((1-rdm)/kaonicFSprob)*sum;//267 361 rdm = ((1-rdm)/kaonicFSprob)*sum;//2670 normalize by the sum of probabilities in the file
362 //now get the line number in the file 362 //now get the line number in the file where the FS particles are stored:
363 G4int n = findStringNumber(rdm, std::m << 363 G4int n = findStringNumber(rdm, probabilities)-1;
364 if ( n < 0 ) return starlist; 364 if ( n < 0 ) return starlist;
365 for(G4int j = 0; j < static_cast<G4int 365 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
366 if(particle_types[n][j] == "pi0"){ 366 if(particle_types[n][j] == "pi0"){
367 Particle *p = new Particle(PiZero, 367 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
368 starlist.push_back(p); 368 starlist.push_back(p);
369 } 369 }
370 else if(particle_types[n][j] == "pi- 370 else if(particle_types[n][j] == "pi-"){
371 Particle *p = new Particle(PiMinus 371 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
372 starlist.push_back(p); 372 starlist.push_back(p);
373 } 373 }
374 else if(particle_types[n][j] == "pi+ 374 else if(particle_types[n][j] == "pi+"){
375 Particle *p = new Particle(PiPlus, 375 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
376 starlist.push_back(p); 376 starlist.push_back(p);
377 } 377 }
378 else if(particle_types[n][j] == "ome 378 else if(particle_types[n][j] == "omega"){
379 Particle *p = new Particle(Omega, 379 Particle *p = new Particle(Omega, mommy, annihilationPosition);
380 starlist.push_back(p); 380 starlist.push_back(p);
381 } 381 }
382 else if(particle_types[n][j] == "eta 382 else if(particle_types[n][j] == "eta"){
383 Particle *p = new Particle(Eta, mo 383 Particle *p = new Particle(Eta, mommy, annihilationPosition);
384 starlist.push_back(p); 384 starlist.push_back(p);
385 } 385 }
386 else if(particle_types[n][j] == "K-" 386 else if(particle_types[n][j] == "K-"){
387 Particle *p = new Particle(KMinus, 387 Particle *p = new Particle(KMinus, mommy, annihilationPosition);
388 starlist.push_back(p); 388 starlist.push_back(p);
389 } 389 }
390 else if(particle_types[n][j] == "K+" 390 else if(particle_types[n][j] == "K+"){
391 Particle *p = new Particle(KPlus, 391 Particle *p = new Particle(KPlus, mommy, annihilationPosition);
392 starlist.push_back(p); 392 starlist.push_back(p);
393 } 393 }
394 else if(particle_types[n][j] == "K0" 394 else if(particle_types[n][j] == "K0"){
395 Particle *p = new Particle(KZero, 395 Particle *p = new Particle(KZero, mommy, annihilationPosition);
396 starlist.push_back(p); 396 starlist.push_back(p);
397 } 397 }
398 else if(particle_types[n][j] == "K0b 398 else if(particle_types[n][j] == "K0b"){
399 Particle *p = new Particle(KZeroBa 399 Particle *p = new Particle(KZeroBar, mommy, annihilationPosition);
400 starlist.push_back(p); 400 starlist.push_back(p);
401 } 401 }
402 else{ 402 else{
403 INCL_ERROR("Some non-existing FS p 403 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
404 for(G4int jj = 0; jj < static_cast 404 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
405 std::cout << "gotcha! " << parti 405 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
406 } 406 }
407 std::cout << "Some non-existing FS 407 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
408 } 408 }
409 } 409 }
410 } 410 }
411 } 411 }
412 else{ //neutronic annihilation 412 else{ //neutronic annihilation
413 INCL_DEBUG("Neutron is the victim" << '\ 413 INCL_DEBUG("Neutron is the victim" << '\n');
414 if(rdm < (1.-kaonicFSprob)){ // pionic/k 414 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
415 INCL_DEBUG("pionic np final state chos 415 INCL_DEBUG("pionic np final state chosen" << '\n');
416 sum = read_file(dataPathnpbar, probabi 416 sum = read_file(dataPathnpbar, probabilities, particle_types);
417 rdm = (rdm/(1.-kaonicFSprob))*sum; //9 417 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.95 normalize by the sum of probabilities in the file
418 //now get the line number in the file 418 //now get the line number in the file where the FS particles are stored:
419 G4int n = findStringNumber(rdm, std::m << 419 G4int n = findStringNumber(rdm, probabilities)-1;
420 if ( n < 0 ) return starlist; 420 if ( n < 0 ) return starlist;
421 for(G4int j = 0; j < static_cast<G4int 421 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
422 if(particle_types[n][j] == "pi0"){ 422 if(particle_types[n][j] == "pi0"){
423 Particle *p = new Particle(PiZero, 423 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
424 starlist.push_back(p); 424 starlist.push_back(p);
425 } 425 }
426 else if(particle_types[n][j] == "pi- 426 else if(particle_types[n][j] == "pi-"){
427 Particle *p = new Particle(PiMinus 427 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
428 starlist.push_back(p); 428 starlist.push_back(p);
429 } 429 }
430 else if(particle_types[n][j] == "pi+ 430 else if(particle_types[n][j] == "pi+"){
431 Particle *p = new Particle(PiPlus, 431 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
432 starlist.push_back(p); 432 starlist.push_back(p);
433 } 433 }
434 else if(particle_types[n][j] == "ome 434 else if(particle_types[n][j] == "omega"){
435 Particle *p = new Particle(Omega, 435 Particle *p = new Particle(Omega, mommy, annihilationPosition);
436 starlist.push_back(p); 436 starlist.push_back(p);
437 } 437 }
438 else if(particle_types[n][j] == "eta 438 else if(particle_types[n][j] == "eta"){
439 Particle *p = new Particle(Eta, mo 439 Particle *p = new Particle(Eta, mommy, annihilationPosition);
440 starlist.push_back(p); 440 starlist.push_back(p);
441 } 441 }
442 else if(particle_types[n][j] == "rho 442 else if(particle_types[n][j] == "rho-"){
443 Particle *p = new Particle(PiMinus 443 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
444 starlist.push_back(p); 444 starlist.push_back(p);
445 Particle *pp = new Particle(PiZero 445 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
446 starlist.push_back(pp); 446 starlist.push_back(pp);
447 } 447 }
448 else if(particle_types[n][j] == "rho 448 else if(particle_types[n][j] == "rho+"){
449 Particle *p = new Particle(PiPlus, 449 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
450 starlist.push_back(p); 450 starlist.push_back(p);
451 Particle *pp = new Particle(PiZero 451 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
452 starlist.push_back(pp); 452 starlist.push_back(pp);
453 } 453 }
454 else if(particle_types[n][j] == "rho 454 else if(particle_types[n][j] == "rho0"){
455 Particle *p = new Particle(PiMinus 455 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
456 starlist.push_back(p); 456 starlist.push_back(p);
457 Particle *pp = new Particle(PiPlus 457 Particle *pp = new Particle(PiPlus, mommy, annihilationPosition);
458 starlist.push_back(pp); 458 starlist.push_back(pp);
459 } 459 }
460 else{ 460 else{
461 INCL_ERROR("Some non-existing FS p 461 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
462 for(G4int jj = 0; jj < static_cast 462 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
463 std::cout << "gotcha! " << parti 463 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
464 } 464 }
465 std::cout << "Some non-existing FS 465 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
466 } 466 }
467 } 467 }
468 } 468 }
469 else{ 469 else{
470 INCL_DEBUG("kaonic np final state chos 470 INCL_DEBUG("kaonic np final state chosen" << '\n');
471 sum = read_file(dataPathnpbark, probab 471 sum = read_file(dataPathnpbark, probabilities, particle_types);
472 rdm = ((1-rdm)/kaonicFSprob)*sum;//383 472 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
473 //now get the line number in the file 473 //now get the line number in the file where the FS particles are stored:
474 G4int n = findStringNumber(rdm, std::m << 474 G4int n = findStringNumber(rdm, probabilities)-1;
475 if ( n < 0 ) return starlist; 475 if ( n < 0 ) return starlist;
476 for(G4int j = 0; j < static_cast<G4int 476 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
477 if(particle_types[n][j] == "pi0"){ 477 if(particle_types[n][j] == "pi0"){
478 Particle *p = new Particle(PiZero, 478 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
479 starlist.push_back(p); 479 starlist.push_back(p);
480 } 480 }
481 else if(particle_types[n][j] == "pi- 481 else if(particle_types[n][j] == "pi-"){
482 Particle *p = new Particle(PiMinus 482 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
483 starlist.push_back(p); 483 starlist.push_back(p);
484 } 484 }
485 else if(particle_types[n][j] == "pi+ 485 else if(particle_types[n][j] == "pi+"){
486 Particle *p = new Particle(PiPlus, 486 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
487 starlist.push_back(p); 487 starlist.push_back(p);
488 } 488 }
489 else if(particle_types[n][j] == "ome 489 else if(particle_types[n][j] == "omega"){
490 Particle *p = new Particle(Omega, 490 Particle *p = new Particle(Omega, mommy, annihilationPosition);
491 starlist.push_back(p); 491 starlist.push_back(p);
492 } 492 }
493 else if(particle_types[n][j] == "eta 493 else if(particle_types[n][j] == "eta"){
494 Particle *p = new Particle(Eta, mo 494 Particle *p = new Particle(Eta, mommy, annihilationPosition);
495 starlist.push_back(p); 495 starlist.push_back(p);
496 } 496 }
497 else if(particle_types[n][j] == "K-" 497 else if(particle_types[n][j] == "K-"){
498 Particle *p = new Particle(KMinus, 498 Particle *p = new Particle(KMinus, mommy, annihilationPosition);
499 starlist.push_back(p); 499 starlist.push_back(p);
500 } 500 }
501 else if(particle_types[n][j] == "K+" 501 else if(particle_types[n][j] == "K+"){
502 Particle *p = new Particle(KPlus, 502 Particle *p = new Particle(KPlus, mommy, annihilationPosition);
503 starlist.push_back(p); 503 starlist.push_back(p);
504 } 504 }
505 else if(particle_types[n][j] == "K0" 505 else if(particle_types[n][j] == "K0"){
506 Particle *p = new Particle(KZero, 506 Particle *p = new Particle(KZero, mommy, annihilationPosition);
507 starlist.push_back(p); 507 starlist.push_back(p);
508 } 508 }
509 else if(particle_types[n][j] == "K0b 509 else if(particle_types[n][j] == "K0b"){
510 Particle *p = new Particle(KZeroBa 510 Particle *p = new Particle(KZeroBar, mommy, annihilationPosition);
511 starlist.push_back(p); 511 starlist.push_back(p);
512 } 512 }
513 else{ 513 else{
514 INCL_ERROR("Some non-existing FS p 514 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
515 for(G4int jj = 0; jj < static_cast 515 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
516 std::cout << "gotcha! " << parti 516 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
517 } 517 }
518 std::cout << "Some non-existing FS 518 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
519 } 519 }
520 } 520 }
521 } 521 }
522 } 522 }
523 523
524 // Correction to the Q-value of the enteri 524 // Correction to the Q-value of the entering particle
525 G4double theCorrection1 = theParticle->get 525 G4double theCorrection1 = theParticle->getEmissionPbarQvalueCorrection(a, z, isProton);
526 G4double theCorrection2 = theParticle->get 526 G4double theCorrection2 = theParticle->getEmissionPbarQvalueCorrection(a, z, !isProton);
527 G4double energyOfMesonStar; 527 G4double energyOfMesonStar;
528 if(isProton == true){ 528 if(isProton == true){
529 energyOfMesonStar = theParticle->getTabl 529 energyOfMesonStar = theParticle->getTableMass() + ParticleTable::getTableMass(a,z,0)
530 -ParticleTable::getTableMass(a-1,z-1,0); 530 -ParticleTable::getTableMass(a-1,z-1,0);
531 } 531 }
532 else{ 532 else{
533 energyOfMesonStar = theParticle->getTabl 533 energyOfMesonStar = theParticle->getTableMass() + ParticleTable::getTableMass(a,z,0)
534 -ParticleTable::getTableMass(a-1,z,0) + 534 -ParticleTable::getTableMass(a-1,z,0) + theCorrection2 - theCorrection1;
535 } 535 }
536 536
537 //compute energies of mesons with a phase- 537 //compute energies of mesons with a phase-space model
538 if(starlist.size() < 2){ 538 if(starlist.size() < 2){
539 INCL_ERROR("should never happen, at leas 539 INCL_ERROR("should never happen, at least 2 final state particles!" << '\n');
540 } 540 }
541 else if(starlist.size() == 2){ 541 else if(starlist.size() == 2){
542 ParticleIter first = starlist.begin(); 542 ParticleIter first = starlist.begin();
543 ParticleIter last = std::next(first, 1); 543 ParticleIter last = std::next(first, 1); //starlist.end() gives an error of segfault, idk why
544 G4double m1 = (*first)->getMass(); 544 G4double m1 = (*first)->getMass();
545 G4double m2 = (*last)->getMass(); 545 G4double m2 = (*last)->getMass();
546 G4double s = energyOfMesonStar*energyOfM 546 G4double s = energyOfMesonStar*energyOfMesonStar;
547 G4double mom1 = std::sqrt(s/4 - (std::po 547 G4double mom1 = std::sqrt(s/4 - (std::pow(m1,2) + std::pow(m2,2))/2 - std::pow(m1,2)*std::pow(m2,2)/s + (std::pow(m1,4) + 2*std::pow(m1*m2,2) + std::pow(m2,4))/(4*s));
548 ThreeVector momentello = Random::normVec 548 ThreeVector momentello = Random::normVector(mom1); //like raffaello :)
549 (*first)->setMomentum(momentello); 549 (*first)->setMomentum(momentello);
550 (*first)->adjustEnergyFromMomentum(); 550 (*first)->adjustEnergyFromMomentum();
551 (*last)->setMomentum(-momentello); 551 (*last)->setMomentum(-momentello);
552 (*last)->adjustEnergyFromMomentum(); 552 (*last)->adjustEnergyFromMomentum();
553 //std::cout << (*first)->getEnergy() << 553 //std::cout << (*first)->getEnergy() << std::endl;
554 } 554 }
555 else{ 555 else{
556 PhaseSpaceGenerator::generate(energyOfMe 556 PhaseSpaceGenerator::generate(energyOfMesonStar, starlist);
557 //ParticleIter first = starlist.begin(); 557 //ParticleIter first = starlist.begin();
558 //std::cout << (*first)->getEnergy() << 558 //std::cout << (*first)->getEnergy() << std::endl;
559 //ParticleIter last = std::next(first, 1 559 //ParticleIter last = std::next(first, 1);
560 //std::cout << (*last)->getEnergy() << s 560 //std::cout << (*last)->getEnergy() << std::endl;
561 } 561 }
562 562
563 return starlist; 563 return starlist;
564 } 564 }
565 565
566 G4bool PbarAtrestEntryChannel::ProtonIsTheVi 566 G4bool PbarAtrestEntryChannel::ProtonIsTheVictim(){
567 if(theNucleus->getAnnihilationType() == PT 567 if(theNucleus->getAnnihilationType() == PType){
568 INCL_DEBUG("isProton" << '\n'); 568 INCL_DEBUG("isProton" << '\n');
569 return true; //proton is annihilated 569 return true; //proton is annihilated
570 } 570 }
571 else if(theNucleus->getAnnihilationType() 571 else if(theNucleus->getAnnihilationType() == NType){
572 INCL_DEBUG("isNeutron" << '\n'); 572 INCL_DEBUG("isNeutron" << '\n');
573 return false; //neutron is annihilated 573 return false; //neutron is annihilated
574 } 574 }
575 else{ 575 else{
576 INCL_ERROR("should never happen, n or p 576 INCL_ERROR("should never happen, n or p is your only choice!" << '\n');
577 G4double rdm3 = Random::shoot(); 577 G4double rdm3 = Random::shoot();
578 if(rdm3 >= 0.){ 578 if(rdm3 >= 0.){
579 // it is set here for test 579 // it is set here for test
580 return false; 580 return false;
581 } 581 }
582 else{ 582 else{
583 return true; 583 return true;
584 } 584 }
585 } 585 }
586 } 586 }
587 587
588 //compute energy lost due to binding with 588 //compute energy lost due to binding with electron shell
589 G4double PbarAtrestEntryChannel::PbarCoulomb 589 G4double PbarAtrestEntryChannel::PbarCoulombicCascadeEnergy(G4int A, G4int Z){
590 G4double N_ann = n_annihilation(A, Z); 590 G4double N_ann = n_annihilation(A, Z);
591 return ParticleTable::getINCLMass(antiProt 591 return ParticleTable::getINCLMass(antiProton)*(A/(A+1.))*(Z*Z/(N_ann*N_ann*2.*137.*137.));
592 592
593 } 593 }
594 /* Coulombic Cascade Energy formula in Boh 594 /* Coulombic Cascade Energy formula in Bohr approximation taken from:
595 Precision spectroscopy of light exotic at 595 Precision spectroscopy of light exotic atoms D. Gotta
596 Progress in Particle and Nuclear Physics 596 Progress in Particle and Nuclear Physics 52 (2004) 133–195
597 This is a crude approximation*/ 597 This is a crude approximation*/
598 598
599 ThreeVector PbarAtrestEntryChannel::getAnnih 599 ThreeVector PbarAtrestEntryChannel::getAnnihilationPosition(){
600 const G4bool isProton = ProtonIsTheVictim( 600 const G4bool isProton = ProtonIsTheVictim();
601 G4int z = theNucleus->getZ(); //was modifi 601 G4int z = theNucleus->getZ(); //was modified in Cascade.cc
602 G4int a = theNucleus->getA(); //was modifi 602 G4int a = theNucleus->getA(); //was modified in Cascade.cc
603 G4double n_ann = n_annihilation(a, z); 603 G4double n_ann = n_annihilation(a, z);
604 a++; //not before the n_ann! 604 a++; //not before the n_ann!
605 605
606 if(isProton == true){z++;} 606 if(isProton == true){z++;}
607 G4double Rpmax = ParticleTable::getMaximum 607 G4double Rpmax = ParticleTable::getMaximumNuclearRadius(Proton, a, z);
608 G4double Rnmax = ParticleTable::getMaximum 608 G4double Rnmax = ParticleTable::getMaximumNuclearRadius(Neutron, a, z);
609 G4double probabilitymax = 0.; //the max va 609 G4double probabilitymax = 0.; //the max value of the probability distribution
610 G4double probability = 0.0; 610 G4double probability = 0.0;
611 G4double radius; 611 G4double radius;
612 612
613 //now we compute the max value of the prob 613 //now we compute the max value of the probability distribution...
614 if(isProton == true){ 614 if(isProton == true){
615 615
616 for(radius = 0.0; radius < Rpmax; radius 616 for(radius = 0.0; radius < Rpmax; radius = radius + 0.001){
617 probability = overlapP(radius, n_ann); 617 probability = overlapP(radius, n_ann);
618 //INCL_WARN("radius, densityP, o 618 //INCL_WARN("radius, densityP, overlapP: " << radius << " " << densityP(radius) << " " << probability << '\n');
619 if(probability > probabilitymax) 619 if(probability > probabilitymax)
620 probabilitymax = probability; //now it 620 probabilitymax = probability; //now it should be the max value of overlapP function
621 } 621 }
622 } 622 }
623 else{ //neutron 623 else{ //neutron
624 624
625 for(radius = 0.0; radius < Rnmax; radius 625 for(radius = 0.0; radius < Rnmax; radius = radius + 0.001){
626 probability = overlapN(radius, n_ann); 626 probability = overlapN(radius, n_ann);
627 //INCL_WARN("radius, densityN, o 627 //INCL_WARN("radius, densityN, overlapN: " << radius << " " << densityN(radius) << " " << probability << '\n');
628 if(probability > probabilitymax) 628 if(probability > probabilitymax)
629 probabilitymax = probability; //now it 629 probabilitymax = probability; //now it should be the max value of overlapP function
630 } 630 }
631 } 631 }
632 632
633 //we know the limits! start rejection algo 633 //we know the limits! start rejection algorithm!
634 G4double x = 0., y = 0.0001, p_for_x = 0.; 634 G4double x = 0., y = 0.0001, p_for_x = 0.;
635 G4double distance = 0.; 635 G4double distance = 0.;
636 if(isProton == true){ 636 if(isProton == true){
637 while(y >= p_for_x){ 637 while(y >= p_for_x){
638 x = Random::shoot() * Rpmax; // create 638 x = Random::shoot() * Rpmax; // create uniformly random r
639 y = Random::shoot() * probabilitymax; 639 y = Random::shoot() * probabilitymax; // create uniformly random prob
640 p_for_x = overlapP(x, n_ann); //probab 640 p_for_x = overlapP(x, n_ann); //probability call for comparison
641 if(y <= p_for_x){ //first cut-off is i 641 if(y <= p_for_x){ //first cut-off is introduced for computational volume reduction
642 distance = x; 642 distance = x;
643 } 643 }
644 } 644 }
645 } 645 }
646 else{ 646 else{
647 while(y >= p_for_x){ 647 while(y >= p_for_x){
648 x = Random::shoot() * Rnmax; // create 648 x = Random::shoot() * Rnmax; // create uniformly random r
649 y = Random::shoot() * probabilitymax; 649 y = Random::shoot() * probabilitymax; // create uniformly random prob
650 p_for_x = overlapN(x, n_ann); //probab 650 p_for_x = overlapN(x, n_ann); //probability call for comparison
651 if(y <= p_for_x){ //first cut-off is i 651 if(y <= p_for_x){ //first cut-off is introduced for computational volume reduction
652 distance = x; 652 distance = x;
653 } 653 }
654 } 654 }
655 } 655 }
656 656
657 //FINAL POSITION VECTOR 657 //FINAL POSITION VECTOR
658 ThreeVector annihilationPosition(0., 0., - 658 ThreeVector annihilationPosition(0., 0., -distance); //3D sphere of distance radius
659 659
660 660
661 return annihilationPosition; 661 return annihilationPosition;
662 } 662 }
663 663
664 664
665 G4double PbarAtrestEntryChannel::n_annihilat 665 G4double PbarAtrestEntryChannel::n_annihilation(G4int A, G4int Z){
666 const G4bool isProton = ProtonIsTheVictim( 666 const G4bool isProton = ProtonIsTheVictim();
667 G4int z = Z; 667 G4int z = Z;
668 G4int a = A; 668 G4int a = A;
669 a++; 669 a++;
670 if(isProton == true){ 670 if(isProton == true){
671 z++; 671 z++;
672 } 672 }
673 INCL_DEBUG("the original Z value is " << z 673 INCL_DEBUG("the original Z value is " << z << '\n');
674 INCL_DEBUG("the original A value is " << a 674 INCL_DEBUG("the original A value is " << a << '\n');
675 G4double n_ann; //annihilation principal q 675 G4double n_ann; //annihilation principal quantum number(interpolation from data H.Poth)
676 if(z <= 1.){ 676 if(z <= 1.){
677 n_ann = 1.; 677 n_ann = 1.;
678 } 678 }
679 else if(z <= 4.){ 679 else if(z <= 4.){
680 n_ann = 2.; 680 n_ann = 2.;
681 } 681 }
682 else if(z <= 11.){ 682 else if(z <= 11.){
683 n_ann = 3.; 683 n_ann = 3.;
684 } 684 }
685 else if(z <= 20.){ 685 else if(z <= 20.){
686 n_ann = 4.; 686 n_ann = 4.;
687 } 687 }
688 else if(z <= 32.){ 688 else if(z <= 32.){
689 n_ann = 5.; 689 n_ann = 5.;
690 } 690 }
691 else if(z <= 46.){ 691 else if(z <= 46.){
692 n_ann = 6.; 692 n_ann = 6.;
693 } 693 }
694 else if(z <= 61.){ 694 else if(z <= 61.){
695 n_ann = 7.; 695 n_ann = 7.;
696 } 696 }
697 else if(z <= 74.){ 697 else if(z <= 74.){
698 n_ann = 8.; 698 n_ann = 8.;
699 } 699 }
700 else if(z <= 84.){ 700 else if(z <= 84.){
701 n_ann = 9.; 701 n_ann = 9.;
702 } 702 }
703 else{ 703 else{
704 n_ann = 10.; 704 n_ann = 10.;
705 } 705 }
706 INCL_DEBUG("The following Pbar will annihi 706 INCL_DEBUG("The following Pbar will annihilate with n = " << n_ann << '\n');
707 707
708 return n_ann; 708 return n_ann;
709 } 709 }
710 710
711 IAvatarList PbarAtrestEntryChannel::bringMes 711 IAvatarList PbarAtrestEntryChannel::bringMesonStar(ParticleList const &pL, Nucleus * const n) {
712 ThreeVector ann_position = getAnnihilation 712 ThreeVector ann_position = getAnnihilationPosition();
713 IAvatarList theAvatarList; 713 IAvatarList theAvatarList;
714 for(ParticleIter p = pL.begin(), e = pL.en 714 for(ParticleIter p = pL.begin(), e = pL.end(); p!=e; ++p){
715 (*p)->setPosition(ann_position); 715 (*p)->setPosition(ann_position);
716 theAvatarList.push_back(new ParticleEntr 716 theAvatarList.push_back(new ParticleEntryAvatar(0.0, n, *p, APAR));
717 } 717 }
718 return theAvatarList; 718 return theAvatarList;
719 } 719 }
720 720
721 void PbarAtrestEntryChannel::fillFinalState( 721 void PbarAtrestEntryChannel::fillFinalState(FinalState *fs) {
722 //const G4bool isProton = ProtonIsTheVicti 722 //const G4bool isProton = ProtonIsTheVictim();
723 //G4int z = theNucleus->getZ(); //was modi 723 //G4int z = theNucleus->getZ(); //was modified in Cascade.cc
724 //G4int a = theNucleus->getA(); //was modi 724 //G4int a = theNucleus->getA(); //was modified in Cascade.cc
725 //a++; //restoration of original A value b 725 //a++; //restoration of original A value before annihilation
726 //if(isProton == true){z++;} //restoration 726 //if(isProton == true){z++;} //restoration of original Z value before annihilation
727 const G4double energyBefore = theParticle- 727 const G4double energyBefore = theParticle->getEnergy();
728 fs->addEnteringParticle(theParticle); 728 fs->addEnteringParticle(theParticle);
729 INCL_DEBUG("Entering particle added " << ' 729 INCL_DEBUG("Entering particle added " << '\n');
730 fs->setTotalEnergyBeforeInteraction(energy 730 fs->setTotalEnergyBeforeInteraction(energyBefore);
731 } 731 }
732 732
733 } 733 }
734 734
735 735
736 736
737 737