Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 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 H 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. 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::PbarAtrestEntryChann 72 :theNucleus(n), theParticle(p) 73 {} 74 75 PbarAtrestEntryChannel::~PbarAtrestEntryChan 76 {} 77 78 //fill probabilities and particle types from 79 G4double PbarAtrestEntryChannel::read_file(s 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::mo 96 } 97 } 98 else std::cout << "ERROR no fread_file " 99 100 return sum_probs; 101 } 102 103 //this function will tell you the FS line nu 104 G4int PbarAtrestEntryChannel::findStringNumb 105 G4int stringNumber = -1; 106 G4double smallestsum = 0.0; 107 G4double biggestsum = yields[0]; 108 //std::cout << "initial input " << rdm << 109 for (G4int i = 0; i < static_cast<G4int>(y 110 if (rdm >= smallestsum && rdm <= bigge 111 //std::cout << smallestsum << " an 112 stringNumber = i+1; 113 } 114 smallestsum += yields[i]; 115 biggestsum += yields[i+1]; 116 } 117 if(stringNumber==-1) stringNumber = static 118 if(stringNumber==-1){ 119 INCL_ERROR("ERROR in findStringNumber (s 120 std::cout << "ERROR in findStringNumber" 121 } 122 return stringNumber; 123 } 124 125 G4double PbarAtrestEntryChannel::fctrl(const 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 G4 133 return std::pow(fctrl(2.0*n),-0.5); 134 } 135 G4double PbarAtrestEntryChannel::r2(const G4 136 G4int Z = theNucleus->getZ(); 137 return std::pow((Z/(14.4*n)), 1.5); 138 } 139 G4double PbarAtrestEntryChannel::r3(G4double 140 G4int Z = theNucleus->getZ(); 141 return std::pow((x)*Z/(n*14.4),(n-1)); //w 142 } 143 G4double PbarAtrestEntryChannel::r4(G4double 144 G4int Z = theNucleus->getZ(); 145 return std::exp((-x*Z)/(n*28.8)); 146 } 147 148 149 G4double PbarAtrestEntryChannel::radial_wave 150 return r1(n)*r2(n)*r3(x,n)*r4(x,n); //Rad 151 } 152 153 /* 154 G4double PbarAtrestEntryChannel::densityP(G4 155 return 0.16800136/(1.0+std::exp((x[0]-par 156 */ 157 G4double PbarAtrestEntryChannel::densityP(G4 158 159 const G4bool isProton = ProtonIsTheVic 160 G4int Z = theNucleus->getZ(); //was mo 161 G4int A = theNucleus->getA(); //was mo 162 A++; //restoration of original A value 163 if(isProton == true){Z++;} //restorati 164 165 if(A > 19) { 166 G4double radius = ParticleTable::get 167 G4double diffuseness = ParticleTable 168 G4double maximumRadius = ParticleTab 169 NuclearDensityFunctions::WoodsSaxon 170 return ((x!=0.) ? rDensityFunction(x 171 } else if(A <= 19 && A > 6) { 172 G4double radius = ParticleTable::get 173 G4double diffuseness = ParticleTable 174 G4double maximumRadius = ParticleTab 175 NuclearDensityFunctions::ModifiedHar 176 return ((x!=0.) ? rDensityFunction(x 177 } else if(A <= 6 && A > 2) { // Gaussi 178 G4double radius = ParticleTable::get 179 G4double maximumRadius = ParticleTab 180 NuclearDensityFunctions::Gaussian rD 181 return ((x!=0.) ? rDensityFunction(x 182 } else if(A == 2 && Z == 1) { // densi 183 NuclearDensityFunctions::ParisR rDen 184 return ((x!=0.) ? rDensityFunction(x 185 } else { 186 INCL_ERROR("No nuclear density funct 187 << A << " Z = " << Z << '\n'); 188 return 0.0; 189 } 190 191 } 192 /* 193 } 194 G4double PbarAtrestEntryChannel::densityN(G4 195 return 0.16800136/(1.0+std::exp((x[0]-par 196 } 197 */ 198 G4double PbarAtrestEntryChannel::densityN(G4 199 200 const G4bool isProton = ProtonIsTheVic 201 G4int Z = theNucleus->getZ(); //was mo 202 G4int A = theNucleus->getA(); //was mo 203 A++; //restoration of original A value 204 if(isProton == true){Z++;} //restorati 205 206 if(A > 19) { 207 G4double radius = ParticleTable::get 208 G4double diffuseness = ParticleTable 209 G4double maximumRadius = ParticleTab 210 NuclearDensityFunctions::WoodsSaxon 211 return ((x!=0.) ? rDensityFunction(x 212 } else if(A <= 19 && A > 6) { 213 G4double radius = ParticleTable::get 214 G4double diffuseness = ParticleTable 215 G4double maximumRadius = ParticleTab 216 NuclearDensityFunctions::ModifiedHar 217 return ((x!=0.) ? rDensityFunction(x 218 } else if(A <= 6 && A > 2) { // Gaussi 219 G4double radius = ParticleTable::get 220 G4double maximumRadius = ParticleTab 221 NuclearDensityFunctions::Gaussian rD 222 return ((x!=0.) ? rDensityFunction(x 223 } else if(A == 2 && Z == 1) { // densi 224 NuclearDensityFunctions::ParisR rDen 225 return ((x!=0.) ? rDensityFunction(x 226 } else { 227 INCL_ERROR("No nuclear density funct 228 << A << " Z = " << Z << '\n'); 229 return 0.0; 230 } 231 232 } 233 234 235 G4double PbarAtrestEntryChannel::overlapP(G4 236 return x*x*r1(n)*r2(n)*r3(x,n)*r4(x,n)*r1 237 } 238 G4double PbarAtrestEntryChannel::overlapN(G4 239 return x*x*r1(n)*r2(n)*r3(x,n)*r4(x,n)*r1 240 } 241 242 243 ParticleList PbarAtrestEntryChannel::makeMes 244 245 // File names 246 #ifdef INCLXX_IN_GEANT4_MODE 247 if(!G4FindDataDir("G4INCLDATA")) { 248 G4ExceptionDescription ed; 249 ed << " Data missing: set environment 250 << " to point to the directory cont 251 << " by the INCL++ model" << G4endl 252 G4Exception("G4INCLDataFile::readDa 253 FatalException, ed); 254 } 255 const G4String& dataPath0{G4FindDataDir( 256 const G4String& dataPathppbar(dataPath0 257 const G4String& dataPathnpbar(dataPath0 258 const G4String& dataPathppbark(dataPath0 259 const G4String& dataPathnpbark(dataPath0 260 #else 261 Config const *theConfig=theNucleus->getS 262 std::string path; 263 if(theConfig) 264 path = theConfig->getINCLXXDataFilePat 265 const std::string& dataPathppbar(path + 266 INCL_DEBUG("Reading https://doi.org/10.1 267 const std::string& dataPathnpbar(path + 268 INCL_DEBUG("Reading https://doi.org/10.1 269 const std::string& dataPathppbark(path + 270 INCL_DEBUG("Reading https://doi.org/10.1 271 const std::string& dataPathnpbark(path + 272 INCL_DEBUG("Reading https://doi.org/10.1 273 #endif 274 /*std::string path = {"/home/zdemid/INCL/i 275 std::string dataPathppbar(path + "/rawppba 276 INCL_DEBUG("Reading https://doi.org/10.101 277 std::string dataPathnpbar(path + "/rawnpba 278 INCL_DEBUG("Reading https://doi.org/10.101 279 std::string dataPathppbark(path + "/rawppb 280 INCL_DEBUG("Reading https://doi.org/10.101 281 std::string dataPathnpbark(path + "/rawnpb 282 INCL_DEBUG("Reading https://doi.org/10.101 283 */ 284 //read probabilities and particle types fr 285 std::vector<G4double> probabilities; //wil 286 std::vector<std::vector<std::string>> part 287 G4double sum; //will contain a sum of prob 288 G4double kaonicFSprob=0.05; //probability 289 290 const G4bool isProton = ProtonIsTheVictim( 291 G4int z = theNucleus->getZ(); //was modifi 292 G4int a = theNucleus->getA(); //was modifi 293 a++; //restoration of original A value bef 294 if(isProton == true){z++;} //restoration o 295 ThreeVector annihilationPosition; 296 ParticleList starlist; 297 ThreeVector mommy; //momentum to be assign 298 299 //LETS GOOOOOOO!!! 300 G4double rdm = Random::shoot(); 301 if(isProton == true){ //protonic annihilat 302 INCL_DEBUG("Proton is the victim" << '\n 303 if(rdm < (1.-kaonicFSprob)){ // pionic F 304 INCL_DEBUG("pionic pp final state chos 305 sum = read_file(dataPathppbar, probabi 306 rdm = (rdm/(1.-kaonicFSprob))*sum; //9 307 //now get the line number in the file 308 G4int n = findStringNumber(rdm, std::m 309 if ( n < 0 ) return starlist; 310 for(G4int j = 0; j < static_cast<G4int 311 if(particle_types[n][j] == "pi0"){ 312 Particle *p = new Particle(PiZero, 313 starlist.push_back(p); 314 } 315 else if(particle_types[n][j] == "pi- 316 Particle *p = new Particle(PiMinus 317 starlist.push_back(p); 318 } 319 else if(particle_types[n][j] == "pi+ 320 Particle *p = new Particle(PiPlus, 321 starlist.push_back(p); 322 } 323 else if(particle_types[n][j] == "ome 324 Particle *p = new Particle(Omega, 325 starlist.push_back(p); 326 } 327 else if(particle_types[n][j] == "eta 328 Particle *p = new Particle(Eta, mo 329 starlist.push_back(p); 330 } 331 else if(particle_types[n][j] == "rho 332 Particle *p = new Particle(PiMinus 333 starlist.push_back(p); 334 Particle *pp = new Particle(PiZero 335 starlist.push_back(pp); 336 } 337 else if(particle_types[n][j] == "rho 338 Particle *p = new Particle(PiPlus, 339 starlist.push_back(p); 340 Particle *pp = new Particle(PiZero 341 starlist.push_back(pp); 342 } 343 else if(particle_types[n][j] == "rho 344 Particle *p = new Particle(PiMinus 345 starlist.push_back(p); 346 Particle *pp = new Particle(PiPlus 347 starlist.push_back(pp); 348 } 349 else{ 350 INCL_ERROR("Some non-existing FS p 351 for(G4int jj = 0; jj < static_cast 352 std::cout << "gotcha! " << parti 353 } 354 std::cout << "Some non-existing FS 355 } 356 } 357 } 358 else{ 359 INCL_DEBUG("kaonic pp final state chos 360 sum = read_file(dataPathppbark, probab 361 rdm = ((1-rdm)/kaonicFSprob)*sum;//267 362 //now get the line number in the file 363 G4int n = findStringNumber(rdm, std::m 364 if ( n < 0 ) return starlist; 365 for(G4int j = 0; j < static_cast<G4int 366 if(particle_types[n][j] == "pi0"){ 367 Particle *p = new Particle(PiZero, 368 starlist.push_back(p); 369 } 370 else if(particle_types[n][j] == "pi- 371 Particle *p = new Particle(PiMinus 372 starlist.push_back(p); 373 } 374 else if(particle_types[n][j] == "pi+ 375 Particle *p = new Particle(PiPlus, 376 starlist.push_back(p); 377 } 378 else if(particle_types[n][j] == "ome 379 Particle *p = new Particle(Omega, 380 starlist.push_back(p); 381 } 382 else if(particle_types[n][j] == "eta 383 Particle *p = new Particle(Eta, mo 384 starlist.push_back(p); 385 } 386 else if(particle_types[n][j] == "K-" 387 Particle *p = new Particle(KMinus, 388 starlist.push_back(p); 389 } 390 else if(particle_types[n][j] == "K+" 391 Particle *p = new Particle(KPlus, 392 starlist.push_back(p); 393 } 394 else if(particle_types[n][j] == "K0" 395 Particle *p = new Particle(KZero, 396 starlist.push_back(p); 397 } 398 else if(particle_types[n][j] == "K0b 399 Particle *p = new Particle(KZeroBa 400 starlist.push_back(p); 401 } 402 else{ 403 INCL_ERROR("Some non-existing FS p 404 for(G4int jj = 0; jj < static_cast 405 std::cout << "gotcha! " << parti 406 } 407 std::cout << "Some non-existing FS 408 } 409 } 410 } 411 } 412 else{ //neutronic annihilation 413 INCL_DEBUG("Neutron is the victim" << '\ 414 if(rdm < (1.-kaonicFSprob)){ // pionic/k 415 INCL_DEBUG("pionic np final state chos 416 sum = read_file(dataPathnpbar, probabi 417 rdm = (rdm/(1.-kaonicFSprob))*sum; //9 418 //now get the line number in the file 419 G4int n = findStringNumber(rdm, std::m 420 if ( n < 0 ) return starlist; 421 for(G4int j = 0; j < static_cast<G4int 422 if(particle_types[n][j] == "pi0"){ 423 Particle *p = new Particle(PiZero, 424 starlist.push_back(p); 425 } 426 else if(particle_types[n][j] == "pi- 427 Particle *p = new Particle(PiMinus 428 starlist.push_back(p); 429 } 430 else if(particle_types[n][j] == "pi+ 431 Particle *p = new Particle(PiPlus, 432 starlist.push_back(p); 433 } 434 else if(particle_types[n][j] == "ome 435 Particle *p = new Particle(Omega, 436 starlist.push_back(p); 437 } 438 else if(particle_types[n][j] == "eta 439 Particle *p = new Particle(Eta, mo 440 starlist.push_back(p); 441 } 442 else if(particle_types[n][j] == "rho 443 Particle *p = new Particle(PiMinus 444 starlist.push_back(p); 445 Particle *pp = new Particle(PiZero 446 starlist.push_back(pp); 447 } 448 else if(particle_types[n][j] == "rho 449 Particle *p = new Particle(PiPlus, 450 starlist.push_back(p); 451 Particle *pp = new Particle(PiZero 452 starlist.push_back(pp); 453 } 454 else if(particle_types[n][j] == "rho 455 Particle *p = new Particle(PiMinus 456 starlist.push_back(p); 457 Particle *pp = new Particle(PiPlus 458 starlist.push_back(pp); 459 } 460 else{ 461 INCL_ERROR("Some non-existing FS p 462 for(G4int jj = 0; jj < static_cast 463 std::cout << "gotcha! " << parti 464 } 465 std::cout << "Some non-existing FS 466 } 467 } 468 } 469 else{ 470 INCL_DEBUG("kaonic np final state chos 471 sum = read_file(dataPathnpbark, probab 472 rdm = ((1-rdm)/kaonicFSprob)*sum;//383 473 //now get the line number in the file 474 G4int n = findStringNumber(rdm, std::m 475 if ( n < 0 ) return starlist; 476 for(G4int j = 0; j < static_cast<G4int 477 if(particle_types[n][j] == "pi0"){ 478 Particle *p = new Particle(PiZero, 479 starlist.push_back(p); 480 } 481 else if(particle_types[n][j] == "pi- 482 Particle *p = new Particle(PiMinus 483 starlist.push_back(p); 484 } 485 else if(particle_types[n][j] == "pi+ 486 Particle *p = new Particle(PiPlus, 487 starlist.push_back(p); 488 } 489 else if(particle_types[n][j] == "ome 490 Particle *p = new Particle(Omega, 491 starlist.push_back(p); 492 } 493 else if(particle_types[n][j] == "eta 494 Particle *p = new Particle(Eta, mo 495 starlist.push_back(p); 496 } 497 else if(particle_types[n][j] == "K-" 498 Particle *p = new Particle(KMinus, 499 starlist.push_back(p); 500 } 501 else if(particle_types[n][j] == "K+" 502 Particle *p = new Particle(KPlus, 503 starlist.push_back(p); 504 } 505 else if(particle_types[n][j] == "K0" 506 Particle *p = new Particle(KZero, 507 starlist.push_back(p); 508 } 509 else if(particle_types[n][j] == "K0b 510 Particle *p = new Particle(KZeroBa 511 starlist.push_back(p); 512 } 513 else{ 514 INCL_ERROR("Some non-existing FS p 515 for(G4int jj = 0; jj < static_cast 516 std::cout << "gotcha! " << parti 517 } 518 std::cout << "Some non-existing FS 519 } 520 } 521 } 522 } 523 524 // Correction to the Q-value of the enteri 525 G4double theCorrection1 = theParticle->get 526 G4double theCorrection2 = theParticle->get 527 G4double energyOfMesonStar; 528 if(isProton == true){ 529 energyOfMesonStar = theParticle->getTabl 530 -ParticleTable::getTableMass(a-1,z-1,0); 531 } 532 else{ 533 energyOfMesonStar = theParticle->getTabl 534 -ParticleTable::getTableMass(a-1,z,0) + 535 } 536 537 //compute energies of mesons with a phase- 538 if(starlist.size() < 2){ 539 INCL_ERROR("should never happen, at leas 540 } 541 else if(starlist.size() == 2){ 542 ParticleIter first = starlist.begin(); 543 ParticleIter last = std::next(first, 1); 544 G4double m1 = (*first)->getMass(); 545 G4double m2 = (*last)->getMass(); 546 G4double s = energyOfMesonStar*energyOfM 547 G4double mom1 = std::sqrt(s/4 - (std::po 548 ThreeVector momentello = Random::normVec 549 (*first)->setMomentum(momentello); 550 (*first)->adjustEnergyFromMomentum(); 551 (*last)->setMomentum(-momentello); 552 (*last)->adjustEnergyFromMomentum(); 553 //std::cout << (*first)->getEnergy() << 554 } 555 else{ 556 PhaseSpaceGenerator::generate(energyOfMe 557 //ParticleIter first = starlist.begin(); 558 //std::cout << (*first)->getEnergy() << 559 //ParticleIter last = std::next(first, 1 560 //std::cout << (*last)->getEnergy() << s 561 } 562 563 return starlist; 564 } 565 566 G4bool PbarAtrestEntryChannel::ProtonIsTheVi 567 if(theNucleus->getAnnihilationType() == PT 568 INCL_DEBUG("isProton" << '\n'); 569 return true; //proton is annihilated 570 } 571 else if(theNucleus->getAnnihilationType() 572 INCL_DEBUG("isNeutron" << '\n'); 573 return false; //neutron is annihilated 574 } 575 else{ 576 INCL_ERROR("should never happen, n or p 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 589 G4double PbarAtrestEntryChannel::PbarCoulomb 590 G4double N_ann = n_annihilation(A, Z); 591 return ParticleTable::getINCLMass(antiProt 592 593 } 594 /* Coulombic Cascade Energy formula in Boh 595 Precision spectroscopy of light exotic at 596 Progress in Particle and Nuclear Physics 597 This is a crude approximation*/ 598 599 ThreeVector PbarAtrestEntryChannel::getAnnih 600 const G4bool isProton = ProtonIsTheVictim( 601 G4int z = theNucleus->getZ(); //was modifi 602 G4int a = theNucleus->getA(); //was modifi 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::getMaximum 608 G4double Rnmax = ParticleTable::getMaximum 609 G4double probabilitymax = 0.; //the max va 610 G4double probability = 0.0; 611 G4double radius; 612 613 //now we compute the max value of the prob 614 if(isProton == true){ 615 616 for(radius = 0.0; radius < Rpmax; radius 617 probability = overlapP(radius, n_ann); 618 //INCL_WARN("radius, densityP, o 619 if(probability > probabilitymax) 620 probabilitymax = probability; //now it 621 } 622 } 623 else{ //neutron 624 625 for(radius = 0.0; radius < Rnmax; radius 626 probability = overlapN(radius, n_ann); 627 //INCL_WARN("radius, densityN, o 628 if(probability > probabilitymax) 629 probabilitymax = probability; //now it 630 } 631 } 632 633 //we know the limits! start rejection algo 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 639 y = Random::shoot() * probabilitymax; 640 p_for_x = overlapP(x, n_ann); //probab 641 if(y <= p_for_x){ //first cut-off is i 642 distance = x; 643 } 644 } 645 } 646 else{ 647 while(y >= p_for_x){ 648 x = Random::shoot() * Rnmax; // create 649 y = Random::shoot() * probabilitymax; 650 p_for_x = overlapN(x, n_ann); //probab 651 if(y <= p_for_x){ //first cut-off is i 652 distance = x; 653 } 654 } 655 } 656 657 //FINAL POSITION VECTOR 658 ThreeVector annihilationPosition(0., 0., - 659 660 661 return annihilationPosition; 662 } 663 664 665 G4double PbarAtrestEntryChannel::n_annihilat 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 674 INCL_DEBUG("the original A value is " << a 675 G4double n_ann; //annihilation principal q 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 annihi 707 708 return n_ann; 709 } 710 711 IAvatarList PbarAtrestEntryChannel::bringMes 712 ThreeVector ann_position = getAnnihilation 713 IAvatarList theAvatarList; 714 for(ParticleIter p = pL.begin(), e = pL.en 715 (*p)->setPosition(ann_position); 716 theAvatarList.push_back(new ParticleEntr 717 } 718 return theAvatarList; 719 } 720 721 void PbarAtrestEntryChannel::fillFinalState( 722 //const G4bool isProton = ProtonIsTheVicti 723 //G4int z = theNucleus->getZ(); //was modi 724 //G4int a = theNucleus->getA(); //was modi 725 //a++; //restoration of original A value b 726 //if(isProton == true){z++;} //restoration 727 const G4double energyBefore = theParticle- 728 fs->addEnteringParticle(theParticle); 729 INCL_DEBUG("Entering particle added " << ' 730 fs->setTotalEnergyBeforeInteraction(energy 731 } 732 733 } 734 735 736 737