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