Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCLNbarAtrestEntryChannel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// INCL++ intra-nuclear cascade model
27// Alain Boudard, CEA-Saclay, France
28// Joseph Cugnon, University of Liege, Belgium
29// Jean-Christophe David, CEA-Saclay, France
30// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31// Sylvie Leray, CEA-Saclay, France
32// Davide Mancusi, CEA-Saclay, France
33//
34#define INCLXX_IN_GEANT4_MODE 1
35
36#include "globals.hh"
37
38/*
39 * G4INCLNbarAtrestEntryChannel.cc
40 *
41 * \date Aug 9, 2024
42 * \author Olivier Lourgo
43 */
44#ifdef INCLXX_IN_GEANT4_MODE
45 #include "G4EnvironmentUtils.hh"
46#endif
47#include "G4INCLNbarAtrestEntryChannel.hh"
48#include "G4INCLRootFinder.hh"
49#include "G4INCLIntersection.hh"
50#include "G4INCLCascade.hh"
51#include <algorithm>
52#include "G4INCLParticle.hh"
53#include "G4INCLKinematicsUtils.hh"
54#include "G4INCLBinaryCollisionAvatar.hh"
55#include "G4INCLRandom.hh"
56#include "G4INCLGlobals.hh"
57#include "G4INCLLogger.hh"
58#include <algorithm>
59#include "G4INCLPhaseSpaceGenerator.hh"
60#include <iostream>
61#include <string>
62#include <sstream>
63#include <vector>
64#include "G4INCLHFB.hh"
65#include "G4INCLParticleEntryAvatar.hh"
66#include "G4INCLNuclearDensityFactory.hh"
67#include "G4INCLNDFWoodsSaxon.hh"
68#include "G4INCLNDFModifiedHarmonicOscillator.hh"
69#include "G4INCLNDFGaussian.hh"
70#include "G4INCLNDFParis.hh"
71#include <string>
72#include <vector>
73#include <iostream>
74#include <fstream>
75#include <sstream>
76
77namespace G4INCL{
78
79 NbarAtrestEntryChannel::NbarAtrestEntryChannel(Nucleus *n, Particle *p)
80 :theNucleus(n), theParticle(p)
81 {}
82 NbarAtrestEntryChannel::~NbarAtrestEntryChannel(){}
83
84 G4double NbarAtrestEntryChannel::read_file(std::string filename, std::vector<G4double>& probabilities, std::vector<std::vector<std::string>>& particle_types){
85 std::ifstream file(filename);
86 G4double sum_probs =0.0;
87 if (file.is_open()){
88 std::string line;
89 while(getline(file,line)){
90 std::istringstream iss(line);
91 G4double prob;
92 iss >> prob;
93 sum_probs += prob;
94 probabilities.push_back(prob);
95 std::vector<std::string> types;
96 std::string type;
97 while (iss >>type){
98 types.push_back(type);
99 }
100 particle_types.push_back(types);
101 }
102 }
103 else std::cout << "ERROR no fread_file " << filename << std::endl;
104
105 return sum_probs;
106 }
107
108
109 G4int NbarAtrestEntryChannel::findStringNumber(G4double rdm, std::vector<G4double> yields){
110 G4int stringNumber =-1;
111 G4double smallestsum =0.0;
112 G4double biggestsum = yields[0];
113 for (G4int i=0; i < static_cast<G4int>(yields.size() -1);i++){
114 if (rdm >= smallestsum && rdm <= biggestsum){
115 stringNumber = i+1;
116
117 }
118 smallestsum += yields[i];
119 biggestsum += yields[i+1];
120 }
121 if (stringNumber==-1) stringNumber = static_cast<G4int>(yields.size());
122 if (stringNumber==-1){
123 INCL_ERROR("ERROR in findStringNumber (stringNumber=-1)");
124 std::cout << "ERROR in findStringNumber" << std::endl;
125 }
126 return stringNumber;
127 }
128
129 G4double NbarAtrestEntryChannel::Pabs(G4double x, G4double value){
130 const G4double r = value; // center of the gaussian
131 const G4double sigma = 1;
132 return std::exp(-std::pow(x-r,2)/(2*sigma*sigma));
133 }
134
135 G4double NbarAtrestEntryChannel::densityP(){ // return the r at which the gaussian of the interaction Probability(Pabs) is centered
136 const G4bool isProton = ProtonIsTheVictim();
137 G4int Z = theNucleus->getZ(); //was modified in Cascade.cc
138 G4int A = theNucleus->getA(); //was modified in Cascade.cc
139 G4double threshold_density = 0.10; //the maximum of the interaction probability is taken at 10% of maximum density
140 //https://doi.org/10.1016/0375-9474(82)90352-9 , Nuclear absorption of stopped antiprotons: Multipion-nucleus interactions, Iljinov, Nazaruk, Chigrinov
141 A++; //restoration of original A value before annihilation
142 if(isProton == true){Z++;} //restoration of original Z value before annihilation
143 if(theNucleus->getAnnihilationType()==DNbarPPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType ||
144 theNucleus->getAnnihilationType()==DNbarNPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType){A++;}
145 if(theNucleus->getAnnihilationType()==DNbarPPbarPType || theNucleus->getAnnihilationType()==DNbarNPbarPType){Z++;}
146
147 if(A > 19) {
148 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
149 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Proton, A, Z);
150 G4double r_10 = diffuseness*std::log((1/threshold_density)-1) + radius; //Radius for a Wood-Saxon
151 return r_10;
152 }else if(A <= 19 && A > 6) {
153 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
154 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Proton, A, Z);
155 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
156 NuclearDensityFunctions::ModifiedHarmonicOscillator rDensityFunction(radius, maximumRadius, diffuseness);
157 //double r_10 = 0.01;
158 //while (rDensityFunction(r_10)/(r_10*r_10) > threshold_density*rDensityFunction(0.01)/(0.01*0.01)) {
159 // r_10 = r_10 + maximumRadius/100. ;
160 //}
161 G4double r_min = 0.01;
162 G4double r_max = maximumRadius;
163 G4double r_10 = (r_min + r_max)/2.;
164 while ((rDensityFunction(r_10)/(r_10*r_10) > 0.11*rDensityFunction(0.01)/(0.01*0.01)) ||
165 (rDensityFunction(r_10)/(r_10*r_10) < 0.09*rDensityFunction(0.01)/(0.01*0.01))) {
166 if (rDensityFunction(r_10)/(r_10*r_10) > 0.11*rDensityFunction(0.01)/(0.01*0.01)) {
167 r_min = r_10;
168 }
169 else {
170 r_max = r_10;
171 }
172 r_10 = (r_min + r_max)/2.;
173 }
174 return r_10;
175 }else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
176 G4double radius = ParticleTable::getRadiusParameter(Proton, A, Z);
177 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Proton, A, Z);
178 NuclearDensityFunctions::Gaussian rDensityFunction(maximumRadius, Math::oneOverSqrtThree * radius);
179 //double r_10=std::sqrt(std::pow(Math::oneOverSqrtThree * radius,2)*std::log(2)); //start when the density is half the maximum
180 //while (rDensityFunction(r_10)/(r_10*r_10) > threshold_density*rDensityFunction(0.01)/(0.01*0.01)) {
181 // r_10 = r_10 + maximumRadius/500. ;
182 //}
183 G4double r_min = 0.01;
184 G4double r_max = maximumRadius;
185 G4double r_10 = (r_min + r_max)/2.;
186 while ((rDensityFunction(r_10)/(r_10*r_10) > 0.11*rDensityFunction(0.01)/(0.01*0.01)) ||
187 (rDensityFunction(r_10)/(r_10*r_10) < 0.09*rDensityFunction(0.01)/(0.01*0.01))) {
188 if (rDensityFunction(r_10)/(r_10*r_10) > 0.11*rDensityFunction(0.01)/(0.01*0.01)) {
189 r_min = r_10;
190 }
191 else {
192 r_max = r_10;
193 }
194 r_10 = (r_min + r_max)/2.;
195 }
196 return r_10;
197 }else {
198 INCL_ERROR("No nuclear density function for target A = "
199 << A << " Z = " << Z << '\n');
200 return 0.0;
201 }
202
203 }
204
205 G4double NbarAtrestEntryChannel::densityN(){
206 const G4bool isProton = ProtonIsTheVictim();
207 G4int Z = theNucleus->getZ(); //was modified in Cascade.cc
208 G4int A = theNucleus->getA(); //was modified in Cascade.cc
209 G4double threshold_density = 0.10; //the maximum of the interaction probability is taken at 10% of maximum density
210 //https://doi.org/10.1016/0375-9474(82)90352-9 , Nuclear absorption of stopped antiprotons: Multipion-nucleus interactions, Iljinov, Nazaruk, Chigrinov
211 A++; //restoration of original A value before annihilation
212 if(isProton == true){Z++;} //restoration of original Z value before annihilation
213 if(theNucleus->getAnnihilationType()==DNbarPPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType ||
214 theNucleus->getAnnihilationType()==DNbarNPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType){A++;}
215 if(theNucleus->getAnnihilationType()==DNbarPPbarPType || theNucleus->getAnnihilationType()==DNbarNPbarPType){Z++;}
216
217 if(A > 19) {
218 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
219 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Neutron, A, Z);
220 G4double r_10 = diffuseness*std::log((1/threshold_density)-1) + radius; //Radius for a Wood-Saxon
221 return r_10;
222 } else if(A <= 19 && A > 6) {
223 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
224 G4double diffuseness = ParticleTable::getSurfaceDiffuseness(Neutron, A, Z);
225 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
226 NuclearDensityFunctions::ModifiedHarmonicOscillator rDensityFunction(radius, maximumRadius, diffuseness);
227 G4double r_10 = 0.01;
228 while (rDensityFunction(r_10)/(r_10*r_10) > threshold_density*rDensityFunction(0.01)/(0.01*0.01)) {
229 r_10 = r_10 + maximumRadius/100. ;
230 }
231 return r_10;
232 } else if(A <= 6 && A > 2) { // Gaussian distribution for light nuclei
233 G4double radius = ParticleTable::getRadiusParameter(Neutron, A, Z);
234 G4double maximumRadius = ParticleTable::getMaximumNuclearRadius(Neutron, A, Z);
235 NuclearDensityFunctions::Gaussian rDensityFunction(maximumRadius, Math::oneOverSqrtThree * radius);
236 G4double r_10=std::sqrt(std::pow(Math::oneOverSqrtThree * radius,2)*std::log(2)); //start when the density is half the maximum
237 while (rDensityFunction(r_10)/(r_10*r_10) > threshold_density*rDensityFunction(0.01)/(0.01*0.01)) {
238 r_10 = r_10 + maximumRadius/500. ;
239 }
240 return r_10;
241 } else {
242 INCL_ERROR("No nuclear density function for target A = "
243 << A << " Z = " << Z << '\n');
244 return 0.0;
245 }
246
247
248 }
249
257 ParticleList NbarAtrestEntryChannel::makeMesonStar() {//This function creates a set of mesons with momenta
258 // File names
259 #ifdef INCLXX_IN_GEANT4_MODE
260 if(!G4FindDataDir("G4INCLDATA")) {
261 G4ExceptionDescription ed;
262 ed << " Data missing: set environment variable G4INCLDATA\n"
263 << " to point to the directory containing data files needed\n"
264 << " by the INCL++ model" << G4endl;
265 G4Exception("G4INCLDataFile::readData()","rawppbarFS.dat, ...",
266 FatalException, ed);
267 }
268 G4String dataPath0{G4FindDataDir("G4INCLDATA")};
269 G4String dataPathnbarp(dataPath0 + "/rawnbarpFS.dat");
270 G4String dataPathnbarn(dataPath0 + "/rawnbarnFS.dat");
271 G4String dataPathnbarnk(dataPath0 + "/rawppbarFSkaonic.dat");
272 G4String dataPathnbarpk(dataPath0 + "/rawnbarpFSkaonic.dat");
273 #else
274 Config const *theConfig=theNucleus->getStore()->getConfig();
275 std::string path;
276 if(theConfig)
277 path = theConfig->getINCLXXDataFilePath();
278 std::string dataPathnbarn(path + "/rawnbarnFS.dat");
279 INCL_DEBUG("Reading nbarn final states" << dataPathnbarn << '\n');
280 std::string dataPathnbarp(path + "/rawnbarpFS.dat");
281 INCL_DEBUG("Reading nbarp final states" << dataPathnbarp << '\n');
282 std::string dataPathnbarnk(path + "/rawppbarFSkaonic.dat");
283 INCL_DEBUG("Reading nbarn kaonic final states" << dataPathnbarnk << '\n');
284 std::string dataPathnbarpk(path + "/rawnbarpFSkaonic.dat");
285 INCL_DEBUG("Reading nbarp kaonic final states" << dataPathnbarpk << '\n');
286 #endif
287 //read probabilities and particle types from file
288 std::vector<G4double> probabilities; //will store each FS yield
289 std::vector<std::vector<std::string>> particle_types; //will store particle names
290
291 G4double sum; //will contain a sum of probabilities of all FS in the file
292 G4double kaonicFSprob=0.05; //probability to kave kaonic FS
293
294 const G4bool isProton = ProtonIsTheVictim();
295 G4int z = theNucleus->getZ(); //was modified in Cascade.cc
296 G4int a = theNucleus->getA(); //was modified in Cascade.cc
297 a++; //restoration of original A value before annihilation
298 if(isProton == true){z++;} //restoration of original Z value before annihilation
299 ThreeVector annihilationPosition;
300 ParticleList starlist;
301 ThreeVector mommy; //momentum to be assigned later
302
303 //LETS GOOOOOOO!!!
304 G4double rdm = Random::shoot();
305 if(isProton == true){ //protonic annihilation
306 INCL_DEBUG("Proton is the victim" << '\n');
307 if(rdm < (1.-kaonicFSprob)){ // pionic FS was chosen
308 INCL_DEBUG("pionic pp final state chosen" << '\n');
309 sum = read_file(dataPathnbarp, probabilities, particle_types);
310 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.88 normalize by the sum of probabilities in the file
311 //now get the line number in the file where the FS particles are stored:
312 G4int n = findStringNumber(rdm, probabilities)-1;
313 if ( n < 0 ) return starlist;
314 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
315 if(particle_types[n][j] == "pi0"){
316 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
317 starlist.push_back(p);
318 }
319 else if(particle_types[n][j] == "pi-"){
320 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
321 starlist.push_back(p);
322 }
323 else if(particle_types[n][j] == "pi+"){
324 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
325 starlist.push_back(p);
326 }
327 else if(particle_types[n][j] == "omega"){
328 Particle *p = new Particle(Omega, mommy, annihilationPosition);
329 starlist.push_back(p);
330 }
331 else if(particle_types[n][j] == "eta"){
332 Particle *p = new Particle(Eta, mommy, annihilationPosition);
333 starlist.push_back(p);
334 }
335 else if(particle_types[n][j] == "rho-"){
336 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
337 starlist.push_back(p);
338 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
339 starlist.push_back(pp);
340 }
341 else if(particle_types[n][j] == "rho+"){
342 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
343 starlist.push_back(p);
344 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
345 starlist.push_back(pp);
346 }
347 else if(particle_types[n][j] == "rho0"){
348 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
349 starlist.push_back(p);
350 Particle *pp = new Particle(PiPlus, mommy, annihilationPosition);
351 starlist.push_back(pp);
352 }
353 else{
354 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
355 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
356 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
357 }
358 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
359 }
360 }
361 }
362 else{
363 INCL_DEBUG("kaonic pp final state chosen" << '\n');
364 sum = read_file(dataPathnbarpk, probabilities, particle_types);
365 rdm = ((1-rdm)/kaonicFSprob)*sum;//2670 normalize by the sum of probabilities in the file
366 //now get the line number in the file where the FS particles are stored:
367 G4int n = findStringNumber(rdm, probabilities)-1;
368 if ( n < 0 ) return starlist;
369 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
370 if(particle_types[n][j] == "pi0"){
371 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
372 starlist.push_back(p);
373 }
374 else if(particle_types[n][j] == "pi-"){
375 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
376 starlist.push_back(p);
377 }
378 else if(particle_types[n][j] == "pi+"){
379 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
380 starlist.push_back(p);
381 }
382 else if(particle_types[n][j] == "omega"){
383 Particle *p = new Particle(Omega, mommy, annihilationPosition);
384 starlist.push_back(p);
385 }
386 else if(particle_types[n][j] == "eta"){
387 Particle *p = new Particle(Eta, mommy, annihilationPosition);
388 starlist.push_back(p);
389 }
390 else if(particle_types[n][j] == "K-"){
391 Particle *p = new Particle(KMinus, mommy, annihilationPosition);
392 starlist.push_back(p);
393 }
394 else if(particle_types[n][j] == "K+"){
395 Particle *p = new Particle(KPlus, mommy, annihilationPosition);
396 starlist.push_back(p);
397 }
398 else if(particle_types[n][j] == "K0"){
399 Particle *p = new Particle(KZero, mommy, annihilationPosition);
400 starlist.push_back(p);
401 }
402 else if(particle_types[n][j] == "K0b"){
403 Particle *p = new Particle(KZeroBar, mommy, annihilationPosition);
404 starlist.push_back(p);
405 }
406 else{
407 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
408 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
409 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
410 }
411 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
412 }
413 }
414 }
415 }
416 else{ //neutronic annihilation
417 INCL_DEBUG("Neutron is the victim" << '\n');
418 if(rdm < (1.-kaonicFSprob)){ // pionic/kaonic choice
419 INCL_DEBUG("pionic np final state chosen" << '\n');
420 sum = read_file(dataPathnbarn, probabilities, particle_types);
421 rdm = (rdm/(1.-kaonicFSprob))*sum; //99.95 normalize by the sum of probabilities in the file
422 //now get the line number in the file where the FS particles are stored:
423 G4int n = findStringNumber(rdm, probabilities)-1;
424 if ( n < 0 ) return starlist;
425 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
426 if(particle_types[n][j] == "pi0"){
427 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
428 starlist.push_back(p);
429 }
430 else if(particle_types[n][j] == "pi-"){
431 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
432 starlist.push_back(p);
433 }
434 else if(particle_types[n][j] == "pi+"){
435 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
436 starlist.push_back(p);
437 }
438 else if(particle_types[n][j] == "omega"){
439 Particle *p = new Particle(Omega, mommy, annihilationPosition);
440 starlist.push_back(p);
441 }
442 else if(particle_types[n][j] == "eta"){
443 Particle *p = new Particle(Eta, mommy, annihilationPosition);
444 starlist.push_back(p);
445 }
446 else if(particle_types[n][j] == "rho-"){
447 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
448 starlist.push_back(p);
449 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
450 starlist.push_back(pp);
451 }
452 else if(particle_types[n][j] == "rho+"){
453 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
454 starlist.push_back(p);
455 Particle *pp = new Particle(PiZero, mommy, annihilationPosition);
456 starlist.push_back(pp);
457 }
458 else if(particle_types[n][j] == "rho0"){
459 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
460 starlist.push_back(p);
461 Particle *pp = new Particle(PiPlus, mommy, annihilationPosition);
462 starlist.push_back(pp);
463 }
464 else{
465 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
466 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
467 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
468 }
469 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
470 }
471 }
472 }
473 else{
474 INCL_DEBUG("kaonic np final state chosen" << '\n');
475 sum = read_file(dataPathnbarnk, probabilities, particle_types);
476 rdm = ((1-rdm)/kaonicFSprob)*sum;//3837 normalize by the sum of probabilities in the file
477 //now get the line number in the file where the FS particles are stored:
478 G4int n = findStringNumber(rdm, probabilities)-1;
479 if ( n < 0 ) return starlist;
480 for(G4int j = 0; j < static_cast<G4int>(particle_types[n].size()); j++){
481 if(particle_types[n][j] == "pi0"){
482 Particle *p = new Particle(PiZero, mommy, annihilationPosition);
483 starlist.push_back(p);
484 }
485 else if(particle_types[n][j] == "pi-"){
486 Particle *p = new Particle(PiMinus, mommy, annihilationPosition);
487 starlist.push_back(p);
488 }
489 else if(particle_types[n][j] == "pi+"){
490 Particle *p = new Particle(PiPlus, mommy, annihilationPosition);
491 starlist.push_back(p);
492 }
493 else if(particle_types[n][j] == "omega"){
494 Particle *p = new Particle(Omega, mommy, annihilationPosition);
495 starlist.push_back(p);
496 }
497 else if(particle_types[n][j] == "eta"){
498 Particle *p = new Particle(Eta, mommy, annihilationPosition);
499 starlist.push_back(p);
500 }
501 else if(particle_types[n][j] == "K-"){
502 Particle *p = new Particle(KMinus, mommy, annihilationPosition);
503 starlist.push_back(p);
504 }
505 else if(particle_types[n][j] == "K+"){
506 Particle *p = new Particle(KPlus, mommy, annihilationPosition);
507 starlist.push_back(p);
508 }
509 else if(particle_types[n][j] == "K0"){
510 Particle *p = new Particle(KZero, mommy, annihilationPosition);
511 starlist.push_back(p);
512 }
513 else if(particle_types[n][j] == "K0b"){
514 Particle *p = new Particle(KZeroBar, mommy, annihilationPosition);
515 starlist.push_back(p);
516 }
517 else{
518 INCL_ERROR("Some non-existing FS particle detected when reading pbar FS files");
519 for(G4int jj = 0; jj < static_cast<G4int>(particle_types[n].size()); jj++){
520 std::cout << "gotcha! " << particle_types[n][jj] << std::endl;
521 }
522 std::cout << "Some non-existing FS particle detected when reading pbar FS files" << std::endl;
523 }
524 }
525 }
526 }
527
528 // Correction to the Q-value of the entering particle
529 G4int stra = theNucleus->getS();
530 G4double energyOfMesonStar;
531 if(theNucleus->isNucleusNucleusCollision()==false){//antiNeutron
532 if(isProton == true){
533 energyOfMesonStar = theParticle->getEnergy() + ParticleTable::getTableMass(a,z,stra)
534 -ParticleTable::getTableMass(a-1,z,stra);
535 }
536 else{
537 energyOfMesonStar = theParticle->getEnergy() + ParticleTable::getTableMass(a,z,stra)
538 -ParticleTable::getTableMass(a-1,z,stra);
539 }
540 } else if(theNucleus->isNucleusNucleusCollision()==true){//antiComposite : job is done in the Antinuclei file
541 return starlist;
542 }
543
544 //compute energies of mesons with a phase-space model
545 if(starlist.size() < 2){
546 INCL_ERROR("should never happen, at least 2 final state particles!" << '\n');
547 }
548 else if(starlist.size() == 2){
549 ParticleIter first = starlist.begin();
550 ParticleIter last = std::next(first, 1); //starlist.end() gives an error of segfault, idk why
551 G4double m1 = (*first)->getMass();
552 G4double m2 = (*last)->getMass();
553 G4double s = energyOfMesonStar*energyOfMesonStar;
554 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));
555 ThreeVector momentello = Random::normVector(mom1); //like raffaello :)
556 (*first)->setMomentum(momentello);
557 (*first)->adjustEnergyFromMomentum();
558 (*last)->setMomentum(-momentello);
559 (*last)->adjustEnergyFromMomentum();
560 //std::cout << (*first)->getEnergy() << std::endl;
561 }
562 else{
563 PhaseSpaceGenerator::generate(energyOfMesonStar, starlist);
564 //ParticleIter first = starlist.begin();
565 //std::cout << (*first)->getEnergy() << std::endl;
566 //ParticleIter last = std::next(first, 1);
567 //std::cout << (*last)->getEnergy() << std::endl;
568 }
569
570 return starlist;
571 }
572
573 G4bool NbarAtrestEntryChannel::ProtonIsTheVictim(){
574 if(theNucleus->getAnnihilationType() == PType || theNucleus->getAnnihilationType() == DNbarPPbarPType || theNucleus->getAnnihilationType() == DNbarPPbarNType ){
575 return true; //a proton is annihilated
576 }
577 else if(theNucleus->getAnnihilationType() == NType || theNucleus->getAnnihilationType() == DNbarNPbarPType || theNucleus->getAnnihilationType() == DNbarNPbarNType){
578 return false; // a neutron is annihilated
579 }
580 else{
581 INCL_ERROR("should never happen, n or p is your only choise" << '\n');
582 G4double rdm3 = Random::shoot();
583 if(rdm3 >= 0.){
584 return false;
585 }
586 else{
587 return true;
588 }
589 }
590 }
591 ThreeVector NbarAtrestEntryChannel::getAnnihilationPosition(){
592 const G4bool isProton = ProtonIsTheVictim();
593 G4int z = theNucleus->getZ(); //was modified in Cascade.cc
594 G4int a = theNucleus->getA(); //was modified in Cascade.cc
595 a++;
596
597 if(isProton == true){z++;}
598 if(theNucleus->getAnnihilationType()==DNbarPPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType ||
599 theNucleus->getAnnihilationType()==DNbarNPbarNType || theNucleus->getAnnihilationType()==DNbarNPbarPType){a++;}
600 if(theNucleus->getAnnihilationType()==DNbarPPbarPType || theNucleus->getAnnihilationType()==DNbarNPbarPType){z++;}
601 G4double Rpmax = ParticleTable::getMaximumNuclearRadius(Proton, a, z);
602 G4double Rnmax = ParticleTable::getMaximumNuclearRadius(Neutron, a, z);
603 G4double probabilitymax = 0.; //the max value of the probability distribution
604 G4double probability = 0.0;
605 G4double radius;
606
607 //now we compute the max value of the probability distribution...
608 if(isProton == true){
609
610 for(radius = 0.0; radius < Rpmax; radius = radius + 0.001){
611 probability = overlapP(radius);
612 //INCL_WARN("radius, densityP, overlapP: " << radius << " " << densityP(radius) << " " << probability << '\n');
613 if(probability > probabilitymax)
614 probabilitymax = probability; //now it should be the max value of overlapP function
615 }
616 }
617 else{ //neutron
618
619 for(radius = 0.0; radius < Rnmax; radius = radius + 0.001){
620 probability = overlapN(radius);
621 //INCL_WARN("radius, densityN, overlapN: " << radius << " " << densityN(radius) << " " << probability << '\n');
622 if(probability > probabilitymax)
623 probabilitymax = probability; //now it should be the max value of overlapP function
624 }
625 }
626
627 //we know the limits! start rejection algorithm!
628 G4double x = 0., y = 0.0001, p_for_x = 0.;
629 G4double distance = 0.;
630 if(isProton == true){
631 while(y >= p_for_x){
632 x = Random::shoot() * Rpmax; // create uniformly random r
633 y = Random::shoot() * probabilitymax; // create uniformly random prob
634 p_for_x = overlapP(x); //probability call for comparison
635 if(y <= p_for_x){ //first cut-off is introduced for computational volume reduction
636 distance = x;
637 }
638 }
639 }
640 else{
641 while(y >= p_for_x){
642 x = Random::shoot() * Rnmax; // create uniformly random r
643 y = Random::shoot() * probabilitymax; // create uniformly random prob
644 p_for_x = overlapN(x); //probability call for comparison
645 if(y <= p_for_x){ //first cut-off is introduced for computational volume reduction
646 distance = x;
647 }
648 }
649 }
650
651 //FINAL POSITION VECTOR
652 //ThreeVector annihilationPosition(0., 0., -distance); //3D sphere of distance radius
653 G4double ctheta = (1.-2.*Random::shoot());
654 G4double stheta = std::sqrt(1.-ctheta*ctheta);
655 G4double phi = Math::twoPi*Random::shoot();
656 ThreeVector annihilationPosition(distance*stheta * std::cos(phi), distance*stheta * std::sin(phi), distance*ctheta); //3D sphere of distance radius
657
658 return annihilationPosition;
659 }
660
661 IAvatarList NbarAtrestEntryChannel::bringMesonStar(ParticleList const &pL, Nucleus * const n) {
662 ThreeVector ann_position = getAnnihilationPosition();
663 IAvatarList theAvatarList;
664 for(ParticleIter p = pL.begin(), e = pL.end(); p!=e; ++p){
665 (*p)->setPosition(ann_position);
666 theAvatarList.push_back(new ParticleEntryAvatar(0.0, n, *p, ANAR));
667 }
668 return theAvatarList;
669 }
670 void NbarAtrestEntryChannel::fillFinalState(FinalState *fs) {
671 //const bool isProton = ProtonIsTheVictim();
672 //int z = theNucleus->getZ(); //was modified in Cascade.cc
673 //int a = theNucleus->getA(); //was modified in Cascade.cc
674 //a++; //restoration of original A value before annihilation
675 //if(isProton == true){z++;} //restoration of original Z value before annihilation
676 const G4double energyBefore = theParticle->getEnergy();
677 fs->addEnteringParticle(theParticle);
678 INCL_DEBUG("Entering particle added " << '\n');
679 fs->setTotalEnergyBeforeInteraction(energyBefore);
680 }
681
682}
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
G4INCLIntersection.hh
Simple class for computing intersections between a straight line and a sphere.
INCL_ERROR
#define INCL_ERROR(x)
Definition G4INCLLogger.hh:225
INCL_DEBUG
#define INCL_DEBUG(x)
Definition G4INCLLogger.hh:240
G4INCLNDFGaussian.hh
Class for Gaussian density.
G4INCLNDFModifiedHarmonicOscillator.hh
Class for modified harmonic oscillator density.
G4INCLNDFParis.hh
NDF* class for the deuteron density according to the Paris potential.
G4INCLNDFWoodsSaxon.hh
Class for Woods-Saxon density.
G4INCLRootFinder.hh
Static root-finder algorithm.
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
G4endl
#define G4endl
Definition G4ios.hh:67
G4INCL::Config::getINCLXXDataFilePath
std::string const & getINCLXXDataFilePath() const
Set the ABLAXX datafile path.
Definition G4INCLConfig.hh:234
G4INCL::FinalState::addEnteringParticle
void addEnteringParticle(Particle *p)
Definition G4INCLFinalState.cc:80
G4INCL::FinalState::setTotalEnergyBeforeInteraction
void setTotalEnergyBeforeInteraction(G4double E)
Definition G4INCLFinalState.hh:71
G4INCL::NbarAtrestEntryChannel::read_file
G4double read_file(std::string filename, std::vector< G4double > &probabilities, std::vector< std::vector< std::string > > &particle_types)
Definition G4INCLNbarAtrestEntryChannel.cc:84
G4INCL::NbarAtrestEntryChannel::bringMesonStar
IAvatarList bringMesonStar(ParticleList const &pL, Nucleus *const n)
Definition G4INCLNbarAtrestEntryChannel.cc:661
G4INCL::NbarAtrestEntryChannel::Pabs
G4double Pabs(G4double x, G4double value)
Definition G4INCLNbarAtrestEntryChannel.cc:129
G4INCL::NbarAtrestEntryChannel::findStringNumber
G4int findStringNumber(G4double rdm, std::vector< G4double > yields)
Definition G4INCLNbarAtrestEntryChannel.cc:109
G4INCL::Math::twoPi
const G4double twoPi
Definition G4INCLGlobals.hh:70
G4INCL::Math::oneOverSqrtThree
const G4double oneOverSqrtThree
Definition G4INCLGlobals.hh:74
G4INCL::ParticleTable::getTableMass
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
Definition G4INCLParticleTable.cc:1492
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::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
G4INCL::Random::shoot
G4double shoot()
Definition G4INCLRandom.cc:93
G4INCL::ParticleIter
ParticleList::const_iterator ParticleIter
Definition G4INCLParticle.hh:72
G4INCL::IAvatarList
UnorderedVector< IAvatar * > IAvatarList
Definition G4INCLIAvatar.hh:108
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