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