Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLPbarAtrestEntryChannel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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