Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
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 711 of file G4INCLPbarAtrestEntryChannel.cc.

711 {
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 }
G4INCL::ParticleIter
ParticleList::const_iterator ParticleIter
Definition G4INCLParticle.hh:72
G4INCL::IAvatarList
UnorderedVector< IAvatar * > IAvatarList
Definition G4INCLIAvatar.hh:107

Referenced by G4INCL::StandardPropagationModel::shootAtrest().

◆ densityN()

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

Definition at line 198 of file G4INCLPbarAtrestEntryChannel.cc.

198 {
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}
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::Particle::getZ
G4int getZ() const
Returns the charge number.
Definition G4INCLParticle.hh:479
G4INCL::Particle::getA
G4int getA() const
Returns the baryon number.
Definition G4INCLParticle.hh:476
G4INCL::Math::oneOverSqrtThree
const G4double oneOverSqrtThree
Definition G4INCLGlobals.hh:73
G4INCL::ParticleTable::getRadiusParameter
G4double getRadiusParameter(const ParticleType t, const G4int A, const G4int Z)
Definition G4INCLParticleTable.cc:1197
G4INCL::ParticleTable::getMaximumNuclearRadius
G4double getMaximumNuclearRadius(const ParticleType t, const G4int A, const G4int Z)
Definition G4INCLParticleTable.cc:1249
G4INCL::ParticleTable::getSurfaceDiffuseness
G4double getSurfaceDiffuseness(const ParticleType t, const G4int A, const G4int Z)
Definition G4INCLParticleTable.cc:1263

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

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 721 of file G4INCLPbarAtrestEntryChannel.cc.

721 {
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 }
INCL_DEBUG
#define INCL_DEBUG(x)
Definition G4INCLLogger.hh:240
G4INCL::Particle::getEnergy
G4double getEnergy() const
Definition G4INCLParticle.hh:910

◆ 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 599 of file G4INCLPbarAtrestEntryChannel.cc.

599 {
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 }
G4INCL::PbarAtrestEntryChannel::overlapP
G4double overlapP(G4double &x, const G4int n)
Definition G4INCLPbarAtrestEntryChannel.cc:235
G4INCL::PbarAtrestEntryChannel::overlapN
G4double overlapN(G4double &x, const G4int n)
Definition G4INCLPbarAtrestEntryChannel.cc:238
G4INCL::Random::shoot
G4double shoot()
Definition G4INCLRandom.cc:93

Referenced by bringMesonStar().

◆ makeMesonStar()

ParticleList G4INCL::PbarAtrestEntryChannel::makeMesonStar ( )

Definition at line 243 of file G4INCLPbarAtrestEntryChannel.cc.

243 {//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 G4String dataPath0{G4FindDataDir("G4INCLDATA")};
256 G4String dataPathppbar(dataPath0 + "/rawppbarFS.dat");
257 G4String dataPathnpbar(dataPath0 + "/rawnpbarFS.dat");
258 G4String dataPathppbark(dataPath0 + "/rawppbarFSkaonic.dat");
259 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 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 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 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 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, 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, 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, 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, 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 }
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::Nucleus::getStore
Store * getStore() const
Definition G4INCLNucleus.hh:361
G4INCL::Particle::getEmissionPbarQvalueCorrection
G4double getEmissionPbarQvalueCorrection(const G4int AParent, const G4int ZParent, const G4bool Victim) const
Definition G4INCLParticle.hh:740
G4INCL::Particle::getTableMass
virtual G4double getTableMass() const
Get the tabulated particle mass.
Definition G4INCLParticle.hh:590
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::Store::getConfig
Config const * getConfig()
Definition G4INCLStore.hh:273
G4INCL::ParticleTable::getTableMass
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
Definition G4INCLParticleTable.cc:1437
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

Referenced by G4INCL::StandardPropagationModel::shootAtrest().

◆ n_annihilation()

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

Definition at line 665 of file G4INCLPbarAtrestEntryChannel.cc.

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

Referenced by getAnnihilationPosition(), and PbarCoulombicCascadeEnergy().

◆ overlapN()

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

Definition at line 238 of file G4INCLPbarAtrestEntryChannel.cc.

238 {
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 }
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 235 of file G4INCLPbarAtrestEntryChannel.cc.

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

Referenced by getAnnihilationPosition().

◆ PbarCoulombicCascadeEnergy()

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

Definition at line 589 of file G4INCLPbarAtrestEntryChannel.cc.

589 {
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 }
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:972

◆ ProtonIsTheVictim()

G4bool G4INCL::PbarAtrestEntryChannel::ProtonIsTheVictim ( )

Definition at line 566 of file G4INCLPbarAtrestEntryChannel.cc.

566 {
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 }
G4INCL::Nucleus::getAnnihilationType
AnnihilationType getAnnihilationType() const
Getter for theAnnihilationType.
Definition G4INCLNucleus.hh:490

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

◆ 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(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: