Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCL::PbarAtrestEntryChannel Class Reference

#include <G4INCLPbarAtrestEntryChannel.hh>

Inheritance diagram for G4INCL::PbarAtrestEntryChannel:

Public Member Functions

 PbarAtrestEntryChannel (Nucleus *n, Particle *p)
virtual ~PbarAtrestEntryChannel ()
void fillFinalState (FinalState *fs)
ParticleList makeMesonStar ()
G4double PbarCoulombicCascadeEnergy (G4int A, G4int Z)
G4double n_annihilation (G4int A, G4int Z)
IAvatarList bringMesonStar (ParticleList const &pL, Nucleus *const n)
G4bool ProtonIsTheVictim ()
ThreeVector getAnnihilationPosition ()
G4double fctrl (const G4double arg1)
G4double r1 (const G4int n)
G4double r2 (const G4int n)
G4double r3 (G4double x, const G4int n)
G4double r4 (G4double x, const G4int n)
G4double radial_wavefunction (G4double x, const G4int n)
G4double densityP (G4double x)
G4double densityN (G4double x)
G4double overlapP (G4double &x, const G4int n)
G4double overlapN (G4double &x, const G4int n)
G4double read_file (std::string filename, std::vector< G4double > &probabilities, std::vector< std::vector< std::string > > &particle_types)
G4int findStringNumber (G4double rdm, std::vector< G4double > yields)
Public Member Functions inherited from G4INCL::IChannel
 IChannel ()
virtual ~IChannel ()
FinalState * getFinalState ()

Detailed Description

Definition at line 59 of file G4INCLPbarAtrestEntryChannel.hh.

Constructor & Destructor Documentation

◆ PbarAtrestEntryChannel()

G4INCL::PbarAtrestEntryChannel::PbarAtrestEntryChannel ( Nucleus * n,
Particle * p )

Definition at line 71 of file G4INCLPbarAtrestEntryChannel.cc.

72 :theNucleus(n), theParticle(p)
73 {}

◆ ~PbarAtrestEntryChannel()

G4INCL::PbarAtrestEntryChannel::~PbarAtrestEntryChannel ( )
virtual

Definition at line 75 of file G4INCLPbarAtrestEntryChannel.cc.

76 {}

Member Function Documentation

◆ bringMesonStar()

IAvatarList G4INCL::PbarAtrestEntryChannel::bringMesonStar ( ParticleList const & pL,
Nucleus *const n )

Definition at line 728 of file G4INCLPbarAtrestEntryChannel.cc.

728 {
729 ThreeVector ann_position = getAnnihilationPosition();
730 IAvatarList theAvatarList;
731 for(ParticleIter p = pL.begin(), e = pL.end(); p!=e; ++p){
732 (*p)->setPosition(ann_position);
733 theAvatarList.push_back(new ParticleEntryAvatar(0.0, n, *p, APAR));
734 }
735 return theAvatarList;
736 }
G4INCL::ParticleIter
ParticleList::const_iterator ParticleIter
Definition G4INCLParticle.hh:72
G4INCL::IAvatarList
UnorderedVector< IAvatar * > IAvatarList
Definition G4INCLIAvatar.hh:108

Referenced by G4INCL::StandardPropagationModel::shootAtrest(), and G4INCL::StandardPropagationModel::shootCompositeAtrest().

◆ densityN()

G4double G4INCL::PbarAtrestEntryChannel::densityN ( G4double x)

Definition at line 201 of file G4INCLPbarAtrestEntryChannel.cc.

201 {
202
203 const G4bool isProton = ProtonIsTheVictim();
204 G4int Z = theNucleus->getZ(); //was modified in Cascade.cc
205 G4int A = theNucleus->getA(); //was modified in Cascade.cc
206 A++; //restoration of original A value before annihilation
207 if(isProton == true){Z++;} //restoration of original Z value before annihilation
208 if(theNucleus->getAnnihilationType()==DNbarPPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType ||
209 theNucleus->getAnnihilationType()==DNbarNPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType){A++;}
210 if(theNucleus->getAnnihilationType()==DNbarPPbarPType || theNucleus->getAnnihilationType()==DNbarPPbarNType){Z++;}
211
212 if(A > 19) {
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::WoodsSaxon rDensityFunction(radius, maximumRadius, diffuseness);
217 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
218 } else if(A <= 19 && A > 6) {
219 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
220 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Neutron, A, Z);
221 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
222 NuclearDensityFunctions::ModifiedHarmonicOscillator rDensityFunction(radius, maximumRadius, diffuseness);
223 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
224 } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
225 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
226 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
227 NuclearDensityFunctions::Gaussian rDensityFunction(maximumRadius, Math::oneOverSqrtThree * radius);
228 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
229 } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
230 NuclearDensityFunctions::ParisR rDensityFunction;
231 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
232 } else {
233 INCL_ERROR("No nuclear density function for target A = "
234 << A << " Z = " << Z << '\n');
235 return 0.0;
236 }
237
238}
INCL_ERROR
#define INCL_ERROR(x)
Definition G4INCLLogger.hh:225
G4double
double G4double
Definition G4Types.hh:83
G4bool
bool G4bool
Definition G4Types.hh:86
G4int
int G4int
Definition G4Types.hh:85
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4INCL::Math::oneOverSqrtThree
const G4double oneOverSqrtThree
Definition G4INCLGlobals.hh:74
G4INCL::ParticleTable::getRadiusParameter
G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z)
Definition G4INCLParticleTable.cc:1230
G4INCL::ParticleTable::getMaximumNuclearRadius
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
Definition G4INCLParticleTable.cc:1282
G4INCL::ParticleTable::getSurfaceDiffuseness
G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z)
Definition G4INCLParticleTable.cc:1296
G4INCL::DNbarNPbarNType
@ DNbarNPbarNType
Definition G4INCLNucleus.hh:65
G4INCL::DNbarNPbarPType
@ DNbarNPbarPType
Definition G4INCLNucleus.hh:65
G4INCL::DNbarPPbarPType
@ DNbarPPbarPType
Definition G4INCLNucleus.hh:65
G4INCL::DNbarPPbarNType
@ DNbarPPbarNType
Definition G4INCLNucleus.hh:65

Referenced by overlapN().

◆ densityP()

G4double G4INCL::PbarAtrestEntryChannel::densityP ( G4double x)

Definition at line 157 of file G4INCLPbarAtrestEntryChannel.cc.

157 {
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 if(theNucleus->getAnnihilationType()==DNbarPPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType ||
165 theNucleus->getAnnihilationType()==DNbarNPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType){A++;}
166 if(theNucleus->getAnnihilationType()==DNbarPPbarPType || theNucleus->getAnnihilationType()==DNbarPPbarNType){Z++;}
167
168 if(A > 19) {
169 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
170 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Proton, A, Z);
171 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
172 NuclearDensityFunctions::WoodsSaxon rDensityFunction(radius, maximumRadius, diffuseness);
173 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
174 } else if(A <= 19 && A > 6) {
175 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
176 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Proton, A, Z);
177 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
178 NuclearDensityFunctions::ModifiedHarmonicOscillator rDensityFunction(radius, maximumRadius, diffuseness);
179 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
180 } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
181 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
182 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
183 NuclearDensityFunctions::Gaussian rDensityFunction(maximumRadius, Math::oneOverSqrtThree * radius);
184 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
185 } else if(A == 2 && Z == 1) { // density from the Paris potential for deuterons
186 NuclearDensityFunctions::ParisR rDensityFunction;
187 return ((x!=0.) ? rDensityFunction(x)/(x*x) : 1.);
188 } else {
189 INCL_ERROR("No nuclear density function for target A = "
190 << A << " Z = " << Z << '\n');
191 return 0.0;
192 }
193
194}

Referenced by overlapP().

◆ fctrl()

G4double G4INCL::PbarAtrestEntryChannel::fctrl ( const G4double arg1)

Definition at line 125 of file G4INCLPbarAtrestEntryChannel.cc.

125 {
126 G4double factorial=1.0;
127 for(G4int k=1; k<=arg1; k++){
128 factorial *= k;
129 }
130 return factorial;
131 }

Referenced by r1().

◆ fillFinalState()

void G4INCL::PbarAtrestEntryChannel::fillFinalState ( FinalState * fs)
virtual

Implements G4INCL::IChannel.

Definition at line 738 of file G4INCLPbarAtrestEntryChannel.cc.

738 {
739 //const G4bool isProton = ProtonIsTheVictim();
740 //G4int z = theNucleus->getZ(); //was modified in Cascade.cc
741 //G4int a = theNucleus->getA(); //was modified in Cascade.cc
742 //a++; //restoration of original A value before annihilation
743 //if(isProton == true){z++;} //restoration of original Z value before annihilation
744 const G4double energyBefore = theParticle->getEnergy();
745 fs->addEnteringParticle(theParticle);
746 INCL_DEBUG("Entering particle added " << '\n');
747 fs->setTotalEnergyBeforeInteraction(energyBefore);
748 }
INCL_DEBUG
#define INCL_DEBUG(x)
Definition G4INCLLogger.hh:240

◆ findStringNumber()

G4int G4INCL::PbarAtrestEntryChannel::findStringNumber ( G4double rdm,
std::vector< G4double > yields )

Definition at line 104 of file G4INCLPbarAtrestEntryChannel.cc.

104 {
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 }

Referenced by makeMesonStar().

◆ getAnnihilationPosition()

ThreeVector G4INCL::PbarAtrestEntryChannel::getAnnihilationPosition ( )

Definition at line 606 of file G4INCLPbarAtrestEntryChannel.cc.

606 {
607 const G4bool isProton = ProtonIsTheVictim();
608 G4int z = theNucleus->getZ(); //was modified in Cascade.cc
609 G4int a = theNucleus->getA(); //was modified in Cascade.cc
610 G4double n_ann = n_annihilation(a, z);
611 a++; //not before the n_ann!
612
613 if(isProton == true){z++;}
614 if(theNucleus->getAnnihilationType()==DNbarPPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType ||
615 theNucleus->getAnnihilationType()==DNbarNPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType){a++;}
616 if(theNucleus->getAnnihilationType()==DNbarPPbarPType || theNucleus->getAnnihilationType()==DNbarPPbarNType){z++;}
617 G4double Rpmax = ParticleTable::getMaximumNuclearRadius(Proton, a, z);
618 G4double Rnmax = ParticleTable::getMaximumNuclearRadius(Neutron, a, z);
619 G4double probabilitymax = 0.; //the max value of the probability distribution
620 G4double probability = 0.0;
621 G4double radius;
622
623 //now we compute the max value of the probability distribution...
624 if(isProton == true){
625
626 for(radius = 0.0; radius < Rpmax; radius = radius + 0.001){
627 probability = overlapP(radius, n_ann);
628 //INCL_WARN("radius, densityP, overlapP: " << radius << " " << densityP(radius) << " " << probability << '\n');
629 if(probability > probabilitymax)
630 probabilitymax = probability; //now it should be the max value of overlapP function
631 }
632 }
633 else{ //neutron
634
635 for(radius = 0.0; radius < Rnmax; radius = radius + 0.001){
636 probability = overlapN(radius, n_ann);
637 //INCL_WARN("radius, densityN, overlapN: " << radius << " " << densityN(radius) << " " << probability << '\n');
638 if(probability > probabilitymax)
639 probabilitymax = probability; //now it should be the max value of overlapP function
640 }
641 }
642
643 //we know the limits! start rejection algorithm!
644 G4double x = 0., y = 0.0001, p_for_x = 0.;
645 G4double distance = 0.;
646 if(isProton == true){
647 while(y >= p_for_x){
648 x = Random::shoot() * Rpmax; // create uniformly random r
649 y = Random::shoot() * probabilitymax; // create uniformly random prob
650 p_for_x = overlapP(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 else{
657 while(y >= p_for_x){
658 x = Random::shoot() * Rnmax; // create uniformly random r
659 y = Random::shoot() * probabilitymax; // create uniformly random prob
660 p_for_x = overlapN(x, n_ann); //probability call for comparison
661 if(y <= p_for_x){ //first cut-off is introduced for computational volume reduction
662 distance = x;
663 }
664 }
665 }
666
667 //FINAL POSITION VECTOR
668 //ThreeVector annihilationPosition(0., 0., -distance); //3D sphere of distance radius
669 G4double ctheta = (1.-2.*Random::shoot());
670 G4double stheta = std::sqrt(1.-ctheta*ctheta);
671 G4double phi = Math::twoPi*Random::shoot();
672 ThreeVector annihilationPosition(distance*stheta * std::cos(phi), distance*stheta * std::sin(phi), distance*ctheta); //3D sphere of distance radius
673
674
675 return annihilationPosition;
676 }
G4INCL::PbarAtrestEntryChannel::overlapP
G4double overlapP(G4double &x, const G4int n)
Definition G4INCLPbarAtrestEntryChannel.cc:241
G4INCL::PbarAtrestEntryChannel::overlapN
G4double overlapN(G4double &x, const G4int n)
Definition G4INCLPbarAtrestEntryChannel.cc:244
G4INCL::Math::twoPi
const G4double twoPi
Definition G4INCLGlobals.hh:70
G4INCL::Random::shoot
G4double shoot()
Definition G4INCLRandom.cc:93

Referenced by bringMesonStar(), and G4INCL::AntinucleiAtrestEntryChannel::makeMesonStar().

◆ makeMesonStar()

ParticleList G4INCL::PbarAtrestEntryChannel::makeMesonStar ( )

Definition at line 249 of file G4INCLPbarAtrestEntryChannel.cc.

249 {//This function creates a set of mesons with momenta
250
251 // File names
252 #ifdef INCLXX_IN_GEANT4_MODE
253 if(!G4FindDataDir("G4INCLDATA")) {
254 G4ExceptionDescription ed;
255 ed << " Data missing: set environment variable G4INCLDATA\n"
256 << " to point to the directory containing data files needed\n"
257 << " by the INCL++ model" << G4endl;
258 G4Exception("G4INCLDataFile::readData()","rawppbarFS.dat, ...",
259 FatalException, ed);
260 }
261 const G4String& dataPath0{G4FindDataDir("G4INCLDATA")};
262 const G4String& dataPathppbar(dataPath0 + "/rawppbarFS.dat");
263 const G4String& dataPathnpbar(dataPath0 + "/rawnpbarFS.dat");
264 const G4String& dataPathppbark(dataPath0 + "/rawppbarFSkaonic.dat");
265 const G4String& dataPathnpbark(dataPath0 + "/rawnpbarFSkaonic.dat");
266 #else
267 Config const *theConfig=theNucleus->getStore()->getConfig();
268 std::string path;
269 if(theConfig)
270 path = theConfig->getINCLXXDataFilePath();
271 const std::string& dataPathppbar(path + "/rawppbarFS.dat");
272 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
273 const std::string& dataPathnpbar(path + "/rawnpbarFS.dat");
274 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
275 const std::string& dataPathppbark(path + "/rawppbarFSkaonic.dat");
276 INCL_DEBUG("Reading https://doi.org/10.1016/j.physrep.2005.03.002 ppbar kaonic final states" << dataPathppbark << '\n');
277 const std::string& dataPathnpbark(path + "/rawnpbarFSkaonic.dat");
278 INCL_DEBUG("Reading https://doi.org/10.1007/BF02818764 and https://link.springer.com/article/10.1007/BF02754930 npbar kaonic final states" << dataPathnpbark << '\n');
279 #endif
280 /*std::string path = {"/home/zdemid/INCL/inclcode/data"};
281 std::string dataPathppbar(path + "/rawppbarFS.dat");
282 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar final states" << dataPathppbar << '\n');
283 std::string dataPathnpbar(path + "/rawnpbarFS.dat");
284 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar final states" << dataPathnpbar << '\n');
285 std::string dataPathppbark(path + "/rawppbarFSkaonic.dat");
286 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N ppbar kaonic final states" << dataPathppbark << '\n');
287 std::string dataPathnpbark(path + "/rawnpbarFSkaonic.dat");
288 INCL_DEBUG("Reading https://doi.org/10.1016/0375-9474(92)90362-N npbar kaonic final states" << dataPathnpbark << '\n');
289 */
290 //read probabilities and particle types from file
291 std::vector<G4double> probabilities; //will store each FS yield
292 std::vector<std::vector<std::string>> particle_types; //will store particle names
293 G4double sum; //will contain a sum of probabilities of all FS in the file
294 G4double kaonicFSprob=0.05; //probability to kave kaonic FS
295
296 const G4bool isProton = ProtonIsTheVictim();
297 G4int z = theNucleus->getZ(); //was modified in Cascade.cc
298 G4int a = theNucleus->getA(); //was modified in Cascade.cc
299 a++; //restoration of original A value before annihilation
300 if(isProton == true){z++;} //restoration of original Z value before annihilation
301 ThreeVector annihilationPosition;
302 ParticleList starlist;
303 ThreeVector mommy; //momentum to be assigned later
304
305 //LETS GOOOOOOO!!!
306 G4double rdm = Random::shoot();
307 if(isProton == true){ //protonic annihilation
308 INCL_DEBUG("Proton is the victim" << '\n');
309 if(rdm < (1.-kaonicFSprob)){ // pionic FS was chosen
310 INCL_DEBUG("pionic pp final state chosen" << '\n');
311 sum = read_file(dataPathppbar, probabilities, particle_types);
312 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.88 normalize by the sum of probabilities in the file
313 //now get the line number in the file where the FS particles are stored:
314 G4int n = findStringNumber(rdm, std::move(probabilities))-1;
315 if ( n < 0 ) return starlist;
316 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
317 if(particle_types[n][j] == "pi0"){
318 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
319 starlist.push_back(p);
320 }
321 else if(particle_types[n][j] == "pi-"){
322 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
323 starlist.push_back(p);
324 }
325 else if(particle_types[n][j] == "pi+"){
326 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
327 starlist.push_back(p);
328 }
329 else if(particle_types[n][j] == "omega"){
330 Particle *p = new Particle(Omega, mommy, annihilationPosition);
331 starlist.push_back(p);
332 }
333 else if(particle_types[n][j] == "eta"){
334 Particle *p = new Particle(Eta, mommy, annihilationPosition);
335 starlist.push_back(p);
336 }
337 else if(particle_types[n][j] == "rho-"){
338 Particle *p = new Particle(PiMinus, 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] == "rho+"){
344 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
345 starlist.push_back(p);
346 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
347 starlist.push_back(pp);
348 }
349 else if(particle_types[n][j] == "rho0"){
350 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
351 starlist.push_back(p);
352 Particle *pp = new Particle(PiPlus, mommy, annihilationPosition);
353 starlist.push_back(pp);
354 }
355 else{
356 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
357 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
358 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
359 }
360 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
361 }
362 }
363 }
364 else{
365 INCL_DEBUG("kaonic pp final state chosen" << '\n');
366 sum = read_file(dataPathppbark, probabilities, particle_types);
367 rdm = ((1-rdm)/kaonicFSprob)*sum;//2670 normalize by the sum of probabilities in the file
368 //now get the line number in the file where the FS particles are stored:
369 G4int n = findStringNumber(rdm, std::move(probabilities))-1;
370 if ( n < 0 ) return starlist;
371 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
372 if(particle_types[n][j] == "pi0"){
373 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
374 starlist.push_back(p);
375 }
376 else if(particle_types[n][j] == "pi-"){
377 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
378 starlist.push_back(p);
379 }
380 else if(particle_types[n][j] == "pi+"){
381 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
382 starlist.push_back(p);
383 }
384 else if(particle_types[n][j] == "omega"){
385 Particle *p = new Particle(Omega, mommy, annihilationPosition);
386 starlist.push_back(p);
387 }
388 else if(particle_types[n][j] == "eta"){
389 Particle *p = new Particle(Eta, mommy, annihilationPosition);
390 starlist.push_back(p);
391 }
392 else if(particle_types[n][j] == "K-"){
393 Particle *p = new Particle(KMinus, mommy, annihilationPosition);
394 starlist.push_back(p);
395 }
396 else if(particle_types[n][j] == "K+"){
397 Particle *p = new Particle(KPlus, mommy, annihilationPosition);
398 starlist.push_back(p);
399 }
400 else if(particle_types[n][j] == "K0"){
401 Particle *p = new Particle(KZero, mommy, annihilationPosition);
402 starlist.push_back(p);
403 }
404 else if(particle_types[n][j] == "K0b"){
405 Particle *p = new Particle(KZeroBar, mommy, annihilationPosition);
406 starlist.push_back(p);
407 }
408 else{
409 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
410 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
411 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
412 }
413 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
414 }
415 }
416 }
417 }
418 else{ //neutronic annihilation
419 INCL_DEBUG("Neutron is the victim" << '\n');
420 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
421 INCL_DEBUG("pionic np final state chosen" << '\n');
422 sum = read_file(dataPathnpbar, probabilities, particle_types);
423 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.95 normalize by the sum of probabilities in the file
424 //now get the line number in the file where the FS particles are stored:
425 G4int n = findStringNumber(rdm, std::move(probabilities))-1;
426 if ( n < 0 ) return starlist;
427 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
428 if(particle_types[n][j] == "pi0"){
429 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
430 starlist.push_back(p);
431 }
432 else if(particle_types[n][j] == "pi-"){
433 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
434 starlist.push_back(p);
435 }
436 else if(particle_types[n][j] == "pi+"){
437 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
438 starlist.push_back(p);
439 }
440 else if(particle_types[n][j] == "omega"){
441 Particle *p = new Particle(Omega, mommy, annihilationPosition);
442 starlist.push_back(p);
443 }
444 else if(particle_types[n][j] == "eta"){
445 Particle *p = new Particle(Eta, mommy, annihilationPosition);
446 starlist.push_back(p);
447 }
448 else if(particle_types[n][j] == "rho-"){
449 Particle *p = new Particle(PiMinus, 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] == "rho+"){
455 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
456 starlist.push_back(p);
457 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
458 starlist.push_back(pp);
459 }
460 else if(particle_types[n][j] == "rho0"){
461 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
462 starlist.push_back(p);
463 Particle *pp = new Particle(PiPlus, mommy, annihilationPosition);
464 starlist.push_back(pp);
465 }
466 else{
467 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
468 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
469 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
470 }
471 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
472 }
473 }
474 }
475 else{
476 INCL_DEBUG("kaonic np final state chosen" << '\n');
477 sum = read_file(dataPathnpbark, probabilities, particle_types);
478 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
479 //now get the line number in the file where the FS particles are stored:
480 G4int n = findStringNumber(rdm, std::move(probabilities))-1;
481 if ( n < 0 ) return starlist;
482 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
483 if(particle_types[n][j] == "pi0"){
484 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
485 starlist.push_back(p);
486 }
487 else if(particle_types[n][j] == "pi-"){
488 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
489 starlist.push_back(p);
490 }
491 else if(particle_types[n][j] == "pi+"){
492 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
493 starlist.push_back(p);
494 }
495 else if(particle_types[n][j] == "omega"){
496 Particle *p = new Particle(Omega, mommy, annihilationPosition);
497 starlist.push_back(p);
498 }
499 else if(particle_types[n][j] == "eta"){
500 Particle *p = new Particle(Eta, mommy, annihilationPosition);
501 starlist.push_back(p);
502 }
503 else if(particle_types[n][j] == "K-"){
504 Particle *p = new Particle(KMinus, mommy, annihilationPosition);
505 starlist.push_back(p);
506 }
507 else if(particle_types[n][j] == "K+"){
508 Particle *p = new Particle(KPlus, mommy, annihilationPosition);
509 starlist.push_back(p);
510 }
511 else if(particle_types[n][j] == "K0"){
512 Particle *p = new Particle(KZero, mommy, annihilationPosition);
513 starlist.push_back(p);
514 }
515 else if(particle_types[n][j] == "K0b"){
516 Particle *p = new Particle(KZeroBar, mommy, annihilationPosition);
517 starlist.push_back(p);
518 }
519 else{
520 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
521 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
522 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
523 }
524 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
525 }
526 }
527 }
528 }
529
530 // Correction to the Q-value of the entering particle
531 G4int stra = theNucleus->getS();
532 G4double energyOfMesonStar;
533 if(theNucleus->isNucleusNucleusCollision()==false){//antiProton
534 if(isProton == true){
535 energyOfMesonStar = theParticle->getEnergy() + ParticleTable::getTableMass(a,z,stra)
536 -ParticleTable::getTableMass(a-1,z-1,stra);
537 }
538 else{
539 energyOfMesonStar = theParticle->getEnergy() + ParticleTable::getTableMass(a,z,stra)
540 -ParticleTable::getTableMass(a-1,z-1,stra);
541 }
542 } else if(theNucleus->isNucleusNucleusCollision()==true){//antiComposite : job is done in the Antinuclei file
543 return starlist;
544 }
545
546 //compute energies of mesons with a phase-space model
547 if(starlist.size() < 2){
548 INCL_ERROR("should never happen, at least 2 final state particles!" << '\n');
549 }
550 else if(starlist.size() == 2){
551 ParticleIter first = starlist.begin();
552 ParticleIter last = std::next(first, 1); //starlist.end() gives an error of segfault, idk why
553 G4double m1 = (*first)->getMass();
554 G4double m2 = (*last)->getMass();
555 G4double s = energyOfMesonStar*energyOfMesonStar;
556 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));
557 ThreeVector momentello = Random::normVector(mom1); //like raffaello :)
558 (*first)->setMomentum(momentello);
559 (*first)->adjustEnergyFromMomentum();
560 (*last)->setMomentum(-momentello);
561 (*last)->adjustEnergyFromMomentum();
562 //std::cout << (*first)->getEnergy() << std::endl;
563 }
564 else{
565 PhaseSpaceGenerator::generate(energyOfMesonStar, starlist);
566 //ParticleIter first = starlist.begin();
567 //std::cout << (*first)->getEnergy() << std::endl;
568 //ParticleIter last = std::next(first, 1);
569 //std::cout << (*last)->getEnergy() << std::endl;
570 }
571
572 return starlist;
573 }
G4FindDataDir
const char * G4FindDataDir(const char *)
Definition G4FindDataDir.cc:94
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition G4Exception.hh:40
G4endl
#define G4endl
Definition G4ios.hh:67
G4INCL::PbarAtrestEntryChannel::read_file
G4double read_file(std::string filename, std::vector< G4double > &probabilities, std::vector< std::vector< std::string > > &particle_types)
Definition G4INCLPbarAtrestEntryChannel.cc:79
G4INCL::PbarAtrestEntryChannel::findStringNumber
G4int findStringNumber(G4double rdm, std::vector< G4double > yields)
Definition G4INCLPbarAtrestEntryChannel.cc:104
G4INCL::ParticleTable::getTableMass
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
Definition G4INCLParticleTable.cc:1492
G4INCL::PhaseSpaceGenerator::generate
void generate(const G4double sqrtS, ParticleList &particles)
Generate an event in the CM system.
Definition G4INCLPhaseSpaceGenerator.cc:94
G4INCL::Random::normVector
ThreeVector normVector(G4double norm=1.)
Definition G4INCLRandom.cc:162
PoPI::ParticleList
std::vector< Base * > ParticleList
Definition PoPI.hpp:186

Referenced by G4INCL::AntinucleiAtrestEntryChannel::makeMesonStar(), G4INCL::StandardPropagationModel::shootAtrest(), and G4INCL::StandardPropagationModel::shootCompositeAtrest().

◆ n_annihilation()

G4double G4INCL::PbarAtrestEntryChannel::n_annihilation ( G4int A,
G4int Z )

Definition at line 679 of file G4INCLPbarAtrestEntryChannel.cc.

679 {
680 const G4bool isProton = ProtonIsTheVictim();
681 G4int z = Z;
682 G4int a = A;
683 a++;
684 if(isProton == true){
685 z++;
686 }
687 if(theNucleus->getAnnihilationType()==DNbarPPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType ||
688 theNucleus->getAnnihilationType()==DNbarNPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType){a++;}
689 if(theNucleus->getAnnihilationType()==DNbarPPbarPType || theNucleus->getAnnihilationType()==DNbarPPbarNType){z++;}
690 INCL_DEBUG("the original Z value is " << z << '\n');
691 INCL_DEBUG("the original A value is " << a << '\n');
692 G4double n_ann; //annihilation principal quantum number(interpolation from data H.Poth)
693 if(z <= 1.){
694 n_ann = 1.;
695 }
696 else if(z <= 4.){
697 n_ann = 2.;
698 }
699 else if(z <= 11.){
700 n_ann = 3.;
701 }
702 else if(z <= 20.){
703 n_ann = 4.;
704 }
705 else if(z <= 32.){
706 n_ann = 5.;
707 }
708 else if(z <= 46.){
709 n_ann = 6.;
710 }
711 else if(z <= 61.){
712 n_ann = 7.;
713 }
714 else if(z <= 74.){
715 n_ann = 8.;
716 }
717 else if(z <= 84.){
718 n_ann = 9.;
719 }
720 else{
721 n_ann = 10.;
722 }
723 INCL_DEBUG("The following Pbar will annihilate with n = " << n_ann << '\n');
724
725 return n_ann;
726 }

Referenced by getAnnihilationPosition(), and PbarCoulombicCascadeEnergy().

◆ overlapN()

G4double G4INCL::PbarAtrestEntryChannel::overlapN ( G4double & x,
const G4int n )

Definition at line 244 of file G4INCLPbarAtrestEntryChannel.cc.

244 {
245 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);
246 }
G4INCL::PbarAtrestEntryChannel::r3
G4double r3(G4double x, const G4int n)
Definition G4INCLPbarAtrestEntryChannel.cc:139
G4INCL::PbarAtrestEntryChannel::r4
G4double r4(G4double x, const G4int n)
Definition G4INCLPbarAtrestEntryChannel.cc:143

Referenced by getAnnihilationPosition().

◆ overlapP()

G4double G4INCL::PbarAtrestEntryChannel::overlapP ( G4double & x,
const G4int n )

Definition at line 241 of file G4INCLPbarAtrestEntryChannel.cc.

241 {
242 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);
243 }

Referenced by getAnnihilationPosition().

◆ PbarCoulombicCascadeEnergy()

G4double G4INCL::PbarAtrestEntryChannel::PbarCoulombicCascadeEnergy ( G4int A,
G4int Z )

Definition at line 596 of file G4INCLPbarAtrestEntryChannel.cc.

596 {
597 G4double N_ann = n_annihilation(A, Z);
598 return ParticleTable::getINCLMass(antiProton)*(A/(A+1.))*(Z*Z/(N_ann*N_ann*2.*137.*137.));
599
600 }
G4INCL::ParticleTable::getINCLMass
G4double getINCLMass(const G4int A, const G4int Z, const G4int S)
Get INCL nuclear mass (in MeV/c^2).
Definition G4INCLParticleTable.cc:1003

◆ ProtonIsTheVictim()

G4bool G4INCL::PbarAtrestEntryChannel::ProtonIsTheVictim ( )

Definition at line 575 of file G4INCLPbarAtrestEntryChannel.cc.

575 {
576 if(theNucleus->getAnnihilationType() == PType || theNucleus->getAnnihilationType() == DNbarNPbarPType || theNucleus->getAnnihilationType() == DNbarPPbarPType){
577 return true; //proton is annihilated
578 }
579 else if(theNucleus->getAnnihilationType() == NType || theNucleus->getAnnihilationType() == DNbarNPbarNType || theNucleus->getAnnihilationType() == DNbarPPbarNType){
580 return false; //neutron is annihilated
581 }
582 else{
583 INCL_ERROR("should never happen, n or p is your only choice!" << '\n');
584 G4double rdm3 = Random::shoot();
585 if(rdm3 >= 0.){
586 // it is set here for test
587 return false;
588 }
589 else{
590 return true;
591 }
592 }
593 }

Referenced by densityN(), densityP(), getAnnihilationPosition(), makeMesonStar(), n_annihilation(), G4INCL::StandardPropagationModel::shootAtrest(), and G4INCL::StandardPropagationModel::shootCompositeAtrest().

◆ r1()

G4double G4INCL::PbarAtrestEntryChannel::r1 ( const G4int n)

Definition at line 132 of file G4INCLPbarAtrestEntryChannel.cc.

132 {
133 return std::pow(fctrl(2.0*n),-0.5);
134 }

Referenced by overlapN(), overlapP(), and radial_wavefunction().

◆ r2()

G4double G4INCL::PbarAtrestEntryChannel::r2 ( const G4int n)

Definition at line 135 of file G4INCLPbarAtrestEntryChannel.cc.

135 {
136 G4int Z = theNucleus->getZ();
137 return std::pow((Z/(14.4*n)), 1.5);
138 }

Referenced by overlapN(), overlapP(), and radial_wavefunction().

◆ r3()

G4double G4INCL::PbarAtrestEntryChannel::r3 ( G4double x,
const G4int n )

Definition at line 139 of file G4INCLPbarAtrestEntryChannel.cc.

139 {
140 G4int Z = theNucleus->getZ();
141 return std::pow((x)*Z/(n*14.4),(n-1)); //why x is vector here?
142 }

Referenced by overlapN(), overlapP(), and radial_wavefunction().

◆ r4()

G4double G4INCL::PbarAtrestEntryChannel::r4 ( G4double x,
const G4int n )

Definition at line 143 of file G4INCLPbarAtrestEntryChannel.cc.

143 {
144 G4int Z = theNucleus->getZ();
145 return std::exp((-x*Z)/(n*28.8));
146 }

Referenced by overlapN(), overlapP(), and radial_wavefunction().

◆ radial_wavefunction()

G4double G4INCL::PbarAtrestEntryChannel::radial_wavefunction ( G4double x,
const G4int n )

Definition at line 149 of file G4INCLPbarAtrestEntryChannel.cc.

149 {
150 return r1(n)*r2(n)*r3(x,n)*r4(x,n); //Radial wave function par=(n, Z)
151 }

◆ read_file()

G4double G4INCL::PbarAtrestEntryChannel::read_file ( std::string filename,
std::vector< G4double > & probabilities,
std::vector< std::vector< std::string > > & particle_types )

Definition at line 79 of file G4INCLPbarAtrestEntryChannel.cc.

79 {
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 }

Referenced by makeMesonStar().

The documentation for this class was generated from the following files: