| 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 #include "G4INCLNNbarToAnnihilationChannel.hh"
39 #include "G4INCLKinematicsUtils.hh"
40 #include "G4INCLBinaryCollisionAvatar.hh"
41 #include "G4INCLRandom.hh"
42 #include "G4INCLGlobals.hh"
43 #include "G4INCLLogger.hh"
44 #include <algorithm>
45 #include "G4INCLPhaseSpaceGenerator.hh"
46
47 namespace G4INCL {
48
49 NNbarToAnnihilationChannel::NNbarToAnnihilat
50 :theNucleus(n), particle1(p1), particle2(p
51 {}
52
53 NNbarToAnnihilationChannel::~NNbarToAnnihila
54
55 //fill probabilities and particle types from d
56 G4double NNbarToAnnihilationChannel::read_file
57 std::ifstream file(filename);
58 G4double sum_probs = 0.0;
59 if (file.is_open()) {
60 std::string line;
61 while (getline(file, line)) {
62 std::istringstream iss(line);
63 G4double prob;
64 iss >> prob;
65 sum_probs += prob;
66 probabilities.push_back(prob);
67 std::vector<std::string> types;
68 std::string type;
69 while (iss >> type) {
70 types.push_back(type);
71 }
72 particle_types.push_back(std::move(t
73 }
74 }
75 else std::cout << "ERROR no fread_file " <<
76 return sum_probs;
77 }
78
79 //this function will tell you the FS line numb
80 G4int NNbarToAnnihilationChannel::findStringNu
81 G4int stringNumber = -1;
82 G4double smallestsum = 0.0;
83 G4double biggestsum = yields[0];
84 //std::cout << "initial input " << rdm <<
85 for (G4int i = 0; i < static_cast<G4int>(y
86 if (rdm >= smallestsum && rdm <= bigge
87 //std::cout << smallestsum << " an
88 stringNumber = i+1;
89 }
90 smallestsum += yields[i];
91 biggestsum += yields[i+1];
92 }
93 if(stringNumber==-1) stringNumber = static
94 if(stringNumber==-1){
95 INCL_ERROR("ERROR in findStringNumber (s
96 std::cout << "ERROR in findStringNumber"
97 }
98 return stringNumber;
99 }
100
101
102 void NNbarToAnnihilationChannel::fillFinalStat
103
104 Particle *nucleon;
105 Particle *antinucleon;
106
107 if(particle1->isNucleon()){
108 nucleon = particle1;
109 antinucleon = particle2;
110 }
111 else{
112 nucleon = particle2;
113 antinucleon = particle1;
114 }
115
116 const G4double plab = 0.001*KinematicsUtil
117 const G4double sqrtS = KinematicsUtils::to
118 G4double rdm = Random::shoot();
119
120 const std::vector<G4double> BFMM6 = {66.09
121 const std::vector<G4double> BFMM1 = {119.0
122 const std::vector<G4double> BFMM471 = {108
123
124 //PPbar annihilation xs
125 const std::vector<G4double> PPbar_pip_pim
126 const std::vector<G4double> PPbar_pip_pim_
127 const std::vector<G4double> PPbar_pip_pim_
128 const std::vector<G4double> PPbar_pip_pim_
129 const std::vector<G4double> PPbar_pip_pim_
130 const std::vector<G4double> PPbar_2pip_2pi
131 const std::vector<G4double> PPbar_2pip_2pi
132 const std::vector<G4double> PPbar_2pip_2pi
133 const std::vector<G4double> PPbar_3pip_3pi
134 const std::vector<G4double> PPbar_3pip_3pi
135 const std::vector<G4double> PPbar_3pip_3pi
136 const std::vector<G4double> PPbar_3pip_3pi
137 const std::vector<G4double> PPbar_4pip_4pi
138 const std::vector<G4double> PPbar_4pip_4pi
139
140 //NPbar annihilation xs
141 const std::vector<G4double> NPbar_pip_2pim
142 const std::vector<G4double> NPbar_pip_2pim
143 const std::vector<G4double> NPbar_pip_2pim
144 const std::vector<G4double> NPbar_pip_2pim
145 const std::vector<G4double> NPbar_pip_pim_
146 const std::vector<G4double> NPbar_pip_pim_
147 const std::vector<G4double> NPbar_2pip_3pi
148 const std::vector<G4double> NPbar_2pip_3pi
149
150
151 // File names
152 #ifdef INCLXX_IN_GEANT4_MODE
153 if(!G4FindDataDir("G4INCLDATA")) {
154 G4ExceptionDescription ed;
155 ed << " Data missing: set environm
156 << " to point to the directory
157 << " by the INCL++ model" << G4
158 G4Exception("G4INCLDataFile::re
159 FatalException, ed);
160 }
161 G4String dataPath0{G4FindDataDir("G4
162 G4String dataPathppbar(dataPath0 + "
163 G4String dataPathnpbar(dataPath0 + "
164 G4String dataPathppbark(dataPath0 +
165 G4String dataPathnpbark(dataPath0 +
166 G4String dataPathpnbar(dataPath0 + "
167 G4String dataPathpnbark(dataPath0 +
168 #else
169 //Config *theConfig = new G4INCL::Co
170 //theConfig->setINCLXXDataFilePath(G4I
171 Config const *theConfig=theNucleus->
172 std::string path;
173 if(theConfig)
174 path = theConfig->getINCLXXDataFil
175 std::string dataPathppbar(path + "/i
176 INCL_DEBUG("Reading https://doi.org/
177 std::string dataPathnpbar(path + "/i
178 INCL_DEBUG("Reading https://doi.org/
179 std::string dataPathppbark(path + "/
180 INCL_DEBUG("Reading https://doi.org/
181 std::string dataPathnpbark(path + "/
182 INCL_DEBUG("Reading https://doi.org/
183 std::string dataPathpnbar(path + "/i
184 std::string dataPathpnbark(path + "/
185 #endif
186 /*std::string path = {"/home/zdemid/IN
187 std::string dataPathppbar(path + "/inf
188 INCL_DEBUG("Reading https://doi.org/10
189 std::string dataPathnpbar(path + "/inf
190 INCL_DEBUG("Reading https://doi.org/10
191 std::string dataPathppbark(path + "/in
192 INCL_DEBUG("Reading https://doi.org/10
193 std::string dataPathnpbark(path + "/in
194 INCL_DEBUG("Reading https://doi.org/10
195 every time we remove lines for which w
196 */
197
198 std::vector<G4double> probabilities; //wil
199 std::vector<std::vector<std::string>> part
200 G4double sum; //will contain a sum of prob
201 const G4double kaonicFSprob=0.05; //probab
202
203
204 ParticleList list;
205 //list.push_back(nucleon);
206 //list.push_back(antinucleon);
207 // NNbar will not be in the list because t
208 const ThreeVector &rcol = nucleon->getPosi
209 const ThreeVector zero;
210
211 //setting types of new particles and pushi
212 if(nucleon->getType()==Neutron && antinucl
213 //std::cout << "npbar"<< std::endl;
214 const G4double totalpnbar = KinematicsUt
215 // xs is same for npbar, but the fs has
216
217 if (rdm * totalpnbar < KinematicsUtils::
218 // First condition
219 Particle* p1 = new Particle(PiPlus,
220 Particle* p2 = new Particle(PiMinus,
221 Particle* p3 = new Particle(PiMinus,
222
223 list.push_back(p1);
224 list.push_back(p2);
225 list.push_back(p3);
226 } else if (rdm * totalpnbar < Kinematics
227 KinematicsUt
228 // Second condition
229 Particle* p1 = new Particle(PiPlus,
230 Particle* p2 = new Particle(PiMinus,
231 Particle* p3 = new Particle(PiZero,
232 Particle* p4 = new Particle(PiZero,
233 Particle* p5 = new Particle(PiMinus,
234
235 list.push_back(p1);
236 list.push_back(p2);
237 list.push_back(p3);
238 list.push_back(p4);
239 list.push_back(p5);
240 } else if (rdm * totalpnbar < Kinematics
241 KinematicsUt
242 KinematicsUt
243 // Third condition
244 Particle* p1 = new Particle(PiPlus,
245 Particle* p2 = new Particle(PiMinus,
246 Particle* p3 = new Particle(PiZero,
247 Particle* p4 = new Particle(PiZero,
248 Particle* p5 = new Particle(PiZero,
249 Particle* p6 = new Particle(PiMinus,
250
251 list.push_back(p1);
252 list.push_back(p2);
253 list.push_back(p3);
254 list.push_back(p4);
255 list.push_back(p5);
256 list.push_back(p6);
257 } else if (rdm * totalpnbar < Kinematics
258 KinematicsUt
259 KinematicsUt
260 KinematicsUt
261 // Fourth condition
262 Particle* p1 = new Particle(PiPlus,
263 Particle* p2 = new Particle(PiMinus,
264 Particle* p3 = new Particle(PiZero,
265 Particle* p4 = new Particle(PiMinus,
266
267 list.push_back(p1);
268 list.push_back(p2);
269 list.push_back(p3);
270 list.push_back(p4);
271 } else if (rdm * totalpnbar < Kinematics
272 KinematicsUt
273 KinematicsUt
274 KinematicsUt
275 KinematicsUt
276 // Fifth condition
277 Particle* p1 = new Particle(PiPlus,
278 Particle* p2 = new Particle(PiMinus,
279 Particle* p3 = new Particle(PiZero,
280 Particle* p4 = new Particle(KMinus,
281 Particle* p5 = new Particle(KZero, z
282
283 list.push_back(p1);
284 list.push_back(p2);
285 list.push_back(p3);
286 list.push_back(p4);
287 list.push_back(p5);
288 } else if (rdm * totalpnbar < Kinematics
289 KinematicsUt
290 KinematicsUt
291 KinematicsUt
292 KinematicsUt
293 KinematicsUt
294 // Sixth condition
295 Particle* p1 = new Particle(PiPlus,
296 Particle* p2 = new Particle(PiMinus,
297 Particle* p3 = new Particle(KMinus,
298 Particle* p4 = new Particle(KZero, z
299
300 list.push_back(p1);
301 list.push_back(p2);
302 list.push_back(p3);
303 list.push_back(p4);
304 } else if (rdm * totalpnbar < Kinematics
305 KinematicsUt
306 KinematicsUt
307 KinematicsUt
308 KinematicsUt
309 KinematicsUt
310 KinematicsUt
311 // Seventh condition
312 Particle* p1 = new Particle(PiPlus,
313 Particle* p2 = new Particle(PiPlus,
314 Particle* p3 = new Particle(PiMinus,
315 Particle* p4 = new Particle(PiMinus,
316 Particle* p5 = new Particle(PiMinus,
317 Particle* p6 = new Particle(PiZero,
318
319 list.push_back(p1);
320 list.push_back(p2);
321 list.push_back(p3);
322 list.push_back(p4);
323 list.push_back(p5);
324 list.push_back(p6);
325 } else if (rdm * totalpnbar < Kinematics
326 KinematicsUt
327 KinematicsUt
328 KinematicsUt
329 KinematicsUt
330 KinematicsUt
331 KinematicsUt
332 KinematicsUt
333 // Eighth condition
334 Particle* p1 = new Particle(PiPlus,
335 Particle* p2 = new Particle(PiPlus,
336 Particle* p3 = new Particle(PiMinus,
337 Particle* p4 = new Particle(PiMinus,
338 Particle* p5 = new Particle(PiMinus,
339
340 list.push_back(p1);
341 list.push_back(p2);
342 list.push_back(p3);
343 list.push_back(p4);
344 list.push_back(p5);
345 } else {
346 // Default condition
347 //std::cout << "default condition pnba
348 if(rdm < (1.-kaonicFSprob)){ // pionic
349 INCL_DEBUG("pionic npbar final s
350 sum = read_file(dataPathnpbar, p
351 rdm = (rdm/(1.-kaonicFSprob))*su
352 //now get the line number in the
353 G4int n = findStringNumber(rdm,
354 for(G4int j = 0; j < static_cast
355 if(particle_types[n][j] == "pi
356 Particle *p = new Particle(P
357 list.push_back(p);
358 }
359 else if(particle_types[n][j] =
360 Particle *p = new Particle(P
361 list.push_back(p);
362 }
363 else if(particle_types[n][j] =
364 Particle *p = new Particle(P
365 list.push_back(p);
366 }
367 else if(particle_types[n][j] =
368 Particle *p = new Particle(O
369 list.push_back(p);
370 }
371 else if(particle_types[n][j] =
372 Particle *p = new Particle(E
373 list.push_back(p);
374 }
375 else{
376 INCL_ERROR("Some non-existin
377 for(G4int jj = 0; jj < stati
378 std::cout << "gotcha! " <<
379 }
380 }
381 }
382 } // end of pionic option
383 else{
384 INCL_DEBUG("kaonic npbar final s
385 sum = read_file(dataPathnpbark,
386 rdm = ((1-rdm)/kaonicFSprob)*sum
387 //now get the line number in the
388 G4int n = findStringNumber(rdm,
389 for(G4int j = 0; j < static_cast
390 if(particle_types[n][j] == "pi
391 Particle *p = new Particle
392 list.push_back(p);
393 }
394 else if(particle_types[n][j] =
395 Particle *p = new Particle
396 list.push_back(p);
397 }
398 else if(particle_types[n][j] =
399 Particle *p = new Particle
400 list.push_back(p);
401 }
402 else if(particle_types[n][j] =
403 Particle *p = new Particle
404 list.push_back(p);
405 }
406 else if(particle_types[n][j] =
407 Particle *p = new Particle
408 list.push_back(p);
409 }
410 else if(particle_types[n][j] =
411 Particle *p = new Particle
412 list.push_back(p);
413 }
414 else if(particle_types[n][j] =
415 Particle *p = new Particle
416 list.push_back(p);
417 }
418 else if(particle_types[n][j] =
419 Particle *p = new Particle
420 list.push_back(p);
421 }
422 else if(particle_types[n][j] =
423 Particle *p = new Particle
424 list.push_back(p);
425 }
426 else{
427 INCL_ERROR("Some non-exist
428 for(G4int jj = 0; jj < sta
429 std::cout << "gotcha! "
430 }
431 }
432 }
433 } // end of kaonic option
434 } // end of default annihilation
435
436 }
437 else if(nucleon->getType()==Proton && anti
438 const G4double totalpnbar = KinematicsUt
439 // xs is same for npbar, but the fs has
440
441 if (rdm * totalpnbar < KinematicsUtils::
442 // First condition
443 Particle* p1 = new Particle(PiPlus,
444 Particle* p2 = new Particle(PiMinus,
445 Particle* p3 = new Particle(PiPlus,
446
447 list.push_back(p1);
448 list.push_back(p2);
449 list.push_back(p3);
450 } else if (rdm * totalpnbar < Kinematics
451 KinematicsUt
452 // Second condition
453 Particle* p1 = new Particle(PiPlus,
454 Particle* p2 = new Particle(PiMinus,
455 Particle* p3 = new Particle(PiZero,
456 Particle* p4 = new Particle(PiZero,
457 Particle* p5 = new Particle(PiPlus,
458
459 list.push_back(p1);
460 list.push_back(p2);
461 list.push_back(p3);
462 list.push_back(p4);
463 list.push_back(p5);
464 } else if (rdm * totalpnbar < Kinematics
465 KinematicsUt
466 KinematicsUt
467 // Third condition
468 Particle* p1 = new Particle(PiPlus,
469 Particle* p2 = new Particle(PiMinus,
470 Particle* p3 = new Particle(PiZero,
471 Particle* p4 = new Particle(PiZero,
472 Particle* p5 = new Particle(PiZero,
473 Particle* p6 = new Particle(PiPlus,
474
475 list.push_back(p1);
476 list.push_back(p2);
477 list.push_back(p3);
478 list.push_back(p4);
479 list.push_back(p5);
480 list.push_back(p6);
481 } else if (rdm * totalpnbar < Kinematics
482 KinematicsUt
483 KinematicsUt
484 KinematicsUt
485 // Fourth condition
486 Particle* p1 = new Particle(PiPlus,
487 Particle* p2 = new Particle(PiMinus,
488 Particle* p3 = new Particle(PiZero,
489 Particle* p4 = new Particle(PiPlus,
490
491 list.push_back(p1);
492 list.push_back(p2);
493 list.push_back(p3);
494 list.push_back(p4);
495 } else if (rdm * totalpnbar < Kinematics
496 KinematicsUt
497 KinematicsUt
498 KinematicsUt
499 KinematicsUt
500 // Fifth condition
501 Particle* p1 = new Particle(PiPlus,
502 Particle* p2 = new Particle(PiMinus,
503 Particle* p3 = new Particle(PiZero,
504 Particle* p4 = new Particle(KPlus, z
505 Particle* p5 = new Particle(KZeroBar
506
507 list.push_back(p1);
508 list.push_back(p2);
509 list.push_back(p3);
510 list.push_back(p4);
511 list.push_back(p5);
512 } else if (rdm * totalpnbar < Kinematics
513 KinematicsUt
514 KinematicsUt
515 KinematicsUt
516 KinematicsUt
517 KinematicsUt
518 // Sixth condition
519 Particle* p1 = new Particle(PiPlus,
520 Particle* p2 = new Particle(PiMinus,
521 Particle* p3 = new Particle(KPlus, z
522 Particle* p4 = new Particle(KZeroBar
523
524 list.push_back(p1);
525 list.push_back(p2);
526 list.push_back(p3);
527 list.push_back(p4);
528 } else if (rdm * totalpnbar < Kinematics
529 KinematicsUt
530 KinematicsUt
531 KinematicsUt
532 KinematicsUt
533 KinematicsUt
534 KinematicsUt
535 // Seventh condition
536 Particle* p1 = new Particle(PiPlus,
537 Particle* p2 = new Particle(PiPlus,
538 Particle* p3 = new Particle(PiMinus,
539 Particle* p4 = new Particle(PiMinus,
540 Particle* p5 = new Particle(PiPlus,
541 Particle* p6 = new Particle(PiZero,
542
543 list.push_back(p1);
544 list.push_back(p2);
545 list.push_back(p3);
546 list.push_back(p4);
547 list.push_back(p5);
548 list.push_back(p6);
549 } else if (rdm * totalpnbar < Kinematics
550 KinematicsUt
551 KinematicsUt
552 KinematicsUt
553 KinematicsUt
554 KinematicsUt
555 KinematicsUt
556 KinematicsUt
557 // Eighth condition
558 Particle* p1 = new Particle(PiPlus,
559 Particle* p2 = new Particle(PiPlus,
560 Particle* p3 = new Particle(PiMinus,
561 Particle* p4 = new Particle(PiMinus,
562 Particle* p5 = new Particle(PiPlus,
563
564 list.push_back(p1);
565 list.push_back(p2);
566 list.push_back(p3);
567 list.push_back(p4);
568 list.push_back(p5);
569 } else {
570 // Default condition
571 if(rdm < (1.-kaonicFSprob)){ // pionic
572 INCL_DEBUG("pionic pnbar final s
573 sum = read_file(dataPathpnbar, p
574 rdm = (rdm/(1.-kaonicFSprob))*su
575 //now get the line number in the
576 G4int n = findStringNumber(rdm,
577 for(G4int j = 0; j < static_cast
578 if(particle_types[n][j] == "pi
579 Particle *p = new Particle(P
580 list.push_back(p);
581 }
582 else if(particle_types[n][j] =
583 Particle *p = new Particle(P
584 list.push_back(p);
585 }
586 else if(particle_types[n][j] =
587 Particle *p = new Particle(P
588 list.push_back(p);
589 }
590 else if(particle_types[n][j] =
591 Particle *p = new Particle(O
592 list.push_back(p);
593 }
594 else if(particle_types[n][j] =
595 Particle *p = new Particle(E
596 list.push_back(p);
597 }
598 else{
599 INCL_ERROR("Some non-existin
600 for(G4int jj = 0; jj < stati
601 std::cout << "gotcha! " <<
602 }
603 }
604 }
605 } // end of pionic option
606 else{
607 INCL_DEBUG("kaonic pnbar final s
608 sum = read_file(dataPathnpbark,
609 rdm = ((1-rdm)/kaonicFSprob)*sum
610 //now get the line number in the
611 G4int n = findStringNumber(rdm,
612 for(G4int j = 0; j < static_cast
613 if(particle_types[n][j] == "pi
614 Particle *p = new Particle
615 list.push_back(p);
616 }
617 else if(particle_types[n][j] =
618 Particle *p = new Particle
619 list.push_back(p);
620 }
621 else if(particle_types[n][j] =
622 Particle *p = new Particle
623 list.push_back(p);
624 }
625 else if(particle_types[n][j] =
626 Particle *p = new Particle
627 list.push_back(p);
628 }
629 else if(particle_types[n][j] =
630 Particle *p = new Particle
631 list.push_back(p);
632 }
633 else if(particle_types[n][j] =
634 Particle *p = new Particle
635 list.push_back(p);
636 }
637 else if(particle_types[n][j] =
638 Particle *p = new Particle
639 list.push_back(p);
640 }
641 else if(particle_types[n][j] =
642 Particle *p = new Particle
643 list.push_back(p);
644 }
645 else if(particle_types[n][j] =
646 Particle *p = new Particle
647 list.push_back(p);
648 }
649 else{
650 INCL_ERROR("Some non-exist
651 for(G4int jj = 0; jj < sta
652 std::cout << "gotcha! "
653 }
654 }
655 } // end of kaonic option
656 } // end of default annihilation
657
658 }
659 else{ //ppbar or nnbar
660 //std::cout << "ppbar or nnbar"<< std::e
661 const G4double totalppbar = KinematicsUt
662 // same for nnbar
663
664 if (rdm * totalppbar < KinematicsUtils::
665 // First condition
666 Particle* p1 = new Particle(PiPlus,
667 Particle* p2 = new Particle(PiMinus,
668
669 list.push_back(p1);
670 list.push_back(p2);
671
672
673 } else if (rdm * totalppbar < Kinematics
674 KinematicsUt
675 // Second condition
676 Particle* p1 = new Particle(PiPlus,
677 Particle* p2 = new Particle(PiMinus,
678 Particle* p3 = new Particle(PiZero,
679
680 list.push_back(p1);
681 list.push_back(p2);
682 list.push_back(p3);
683 } else if (rdm * totalppbar < Kinematics
684 KinematicsUt
685 KinematicsUt
686 // Third condition
687 Particle* p1 = new Particle(PiPlus,
688 Particle* p2 = new Particle(PiMinus,
689 Particle* p3 = new Particle(Omega, z
690
691 list.push_back(p1);
692 list.push_back(p2);
693 list.push_back(p3);
694 } else if (rdm * totalppbar < Kinematics
695 KinematicsUt
696 KinematicsUt
697 KinematicsUt
698 // Fourth condition
699 Particle* p1 = new Particle(PiPlus,
700 Particle* p2 = new Particle(PiMinus,
701 Particle* p3 = new Particle(KPlus, z
702 Particle* p4 = new Particle(KMinus,
703
704 list.push_back(p1);
705 list.push_back(p2);
706 list.push_back(p3);
707 list.push_back(p4);
708 } else if (rdm * totalppbar < Kinematics
709 KinematicsUt
710 KinematicsUt
711 KinematicsUt
712 KinematicsUt
713 // Fifth condition
714 Particle* p1 = new Particle(PiPlus,
715 Particle* p2 = new Particle(PiMinus,
716 Particle* p3 = new Particle(PiZero,
717 Particle* p4 = new Particle(KPlus, z
718 Particle* p5 = new Particle(KMinus,
719
720 list.push_back(p1);
721 list.push_back(p2);
722 list.push_back(p3);
723 list.push_back(p4);
724 list.push_back(p5);
725 } else if (rdm * totalppbar < Kinematics
726 KinematicsUt
727 KinematicsUt
728 KinematicsUt
729 KinematicsUt
730 KinematicsUt
731 // Sixth condition
732 Particle* p1 = new Particle(PiPlus,
733 Particle* p2 = new Particle(PiMinus,
734 Particle* p3 = new Particle(PiPlus,
735 Particle* p4 = new Particle(PiMinus,
736
737 list.push_back(p1);
738 list.push_back(p2);
739 list.push_back(p3);
740 list.push_back(p4);
741 } else if (rdm * totalppbar < Kinematics
742 KinematicsUt
743 KinematicsUt
744 KinematicsUt
745 KinematicsUt
746 KinematicsUt
747 KinematicsUt
748 // Seventh condition
749 Particle* p1 = new Particle(PiPlus,
750 Particle* p2 = new Particle(PiMinus,
751 Particle* p3 = new Particle(PiPlus,
752 Particle* p4 = new Particle(PiMinus,
753 Particle* p5 = new Particle(PiZero,
754
755 list.push_back(p1);
756 list.push_back(p2);
757 list.push_back(p3);
758 list.push_back(p4);
759 list.push_back(p5);
760 } else if (rdm * totalppbar < Kinematics
761 KinematicsUtils::c
762 KinematicsUtils::c
763 KinematicsUtils::c
764 KinematicsUtils::c
765 KinematicsUtils::c
766 KinematicsUtils::c
767 KinematicsUtils::c
768 // Eighth condition
769 Particle* p1 = new Particle(PiPlus,
770 Particle* p2 = new Particle(PiMinus,
771 Particle* p3 = new Particle(PiPlus,
772 Particle* p4 = new Particle(PiMinus,
773 Particle* p5 = new Particle(PiZero,
774 Particle* p6 = new Particle(PiZero,
775 Particle* p7 = new Particle(PiZero,
776
777 list.push_back(p1);
778 list.push_back(p2);
779 list.push_back(p3);
780 list.push_back(p4);
781 list.push_back(p5);
782 list.push_back(p6);
783 list.push_back(p7);
784 } else if (rdm * totalppbar < Kinematics
785 KinematicsUt
786 KinematicsUt
787 KinematicsUt
788 KinematicsUt
789 KinematicsUt
790 KinematicsUt
791 KinematicsUt
792 KinematicsUt
793 // Ninth condition
794 Particle* p1 = new Particle(PiPlus,
795 Particle* p2 = new Particle(PiMinus,
796 Particle* p3 = new Particle(PiPlus,
797 Particle* p4 = new Particle(PiMinus,
798 Particle* p5 = new Particle(PiPlus,
799 Particle* p6 = new Particle(PiMinus,
800
801 list.push_back(p1);
802 list.push_back(p2);
803 list.push_back(p3);
804 list.push_back(p4);
805 list.push_back(p5);
806 list.push_back(p6);
807 } else if (rdm * totalppbar < Kinematics
808 KinematicsUt
809 KinematicsUt
810 KinematicsUt
811 KinematicsUt
812 KinematicsUt
813 KinematicsUt
814 KinematicsUt
815 KinematicsUt
816 KinematicsUt
817 // Tenth condition
818 Particle* p1 = new Particle(PiPlus,
819 Particle* p2 = new Particle(PiMinus,
820 Particle* p3 = new Particle(PiPlus,
821 Particle* p4 = new Particle(PiMinus,
822 Particle* p5 = new Particle(PiPlus,
823 Particle* p6 = new Particle(PiMinus,
824 Particle* p7 = new Particle(PiZero,
825
826 list.push_back(p1);
827 list.push_back(p2);
828 list.push_back(p3);
829 list.push_back(p4);
830 list.push_back(p5);
831 list.push_back(p6);
832 list.push_back(p7);
833 } else if (rdm * totalppbar < Kinematics
834 KinematicsUt
835 KinematicsUt
836 KinematicsUt
837 KinematicsUt
838 KinematicsUt
839 KinematicsUt
840 KinematicsUt
841 KinematicsUt
842 KinematicsUt
843 KinematicsUt
844 // Eleventh condition
845 Particle* p1 = new Particle(PiPlus,
846 Particle* p2 = new Particle(PiMinus,
847 Particle* p3 = new Particle(PiPlus,
848 Particle* p4 = new Particle(PiMinus,
849 Particle* p5 = new Particle(PiPlus,
850 Particle* p6 = new Particle(PiMinus,
851 Particle* p7 = new Particle(PiZero,
852 Particle* p8 = new Particle(PiZero,
853
854 list.push_back(p1);
855 list.push_back(p2);
856 list.push_back(p3);
857 list.push_back(p4);
858 list.push_back(p5);
859 list.push_back(p6);
860 list.push_back(p7);
861 list.push_back(p8);
862 } else if (rdm * totalppbar < Kinematics
863 KinematicsUt
864 KinematicsUt
865 KinematicsUt
866 KinematicsUt
867 KinematicsUt
868 KinematicsUt
869 KinematicsUt
870 KinematicsUt
871 KinematicsUt
872 KinematicsUt
873 KinematicsUt
874 // Twelfth condition
875 Particle* p1 = new Particle(PiPlus,
876 Particle* p2 = new Particle(PiMinus,
877 Particle* p3 = new Particle(PiPlus,
878 Particle* p4 = new Particle(PiMinus,
879 Particle* p5 = new Particle(PiPlus,
880 Particle* p6 = new Particle(PiMinus,
881 Particle* p7 = new Particle(PiZero,
882 Particle* p8 = new Particle(PiZero,
883 Particle* p9 = new Particle(PiZero,
884
885 list.push_back(p1);
886 list.push_back(p2);
887 list.push_back(p3);
888 list.push_back(p4);
889 list.push_back(p5);
890 list.push_back(p6);
891 list.push_back(p7);
892 list.push_back(p8);
893 list.push_back(p9);
894 } else if (rdm * totalppbar < Kinematics
895 KinematicsUt
896 KinematicsUt
897 KinematicsUt
898 KinematicsUt
899 KinematicsUt
900 KinematicsUt
901 KinematicsUt
902 KinematicsUt
903 KinematicsUt
904 KinematicsUt
905 KinematicsUt
906 KinematicsUt
907 // Thirteenth condition
908 Particle* p1 = new Particle(PiPlus,
909 Particle* p2 = new Particle(PiMinus,
910 Particle* p3 = new Particle(PiPlus,
911 Particle* p4 = new Particle(PiMinus,
912 Particle* p5 = new Particle(PiPlus,
913 Particle* p6 = new Particle(PiMinus,
914 Particle* p7 = new Particle(PiPlus,
915 Particle* p8 = new Particle(PiMinus,
916
917 list.push_back(p1);
918 list.push_back(p2);
919 list.push_back(p3);
920 list.push_back(p4);
921 list.push_back(p5);
922 list.push_back(p6);
923 list.push_back(p7);
924 list.push_back(p8);
925 } else if (rdm * totalppbar < Kinematics
926 KinematicsUt
927 KinematicsUt
928 KinematicsUt
929 KinematicsUt
930 KinematicsUt
931 KinematicsUt
932 KinematicsUt
933 KinematicsUt
934 KinematicsUt
935 KinematicsUt
936 KinematicsUt
937 KinematicsUt
938 KinematicsUt
939 // Fourteenth condition
940 Particle* p1 = new Particle(PiPlus,
941 Particle* p2 = new Particle(PiMinus,
942 Particle* p3 = new Particle(PiPlus,
943 Particle* p4 = new Particle(PiMinus,
944 Particle* p5 = new Particle(PiPlus,
945 Particle* p6 = new Particle(PiMinus,
946 Particle* p7 = new Particle(PiPlus,
947 Particle* p8 = new Particle(PiMinus,
948 Particle* p9 = new Particle(PiZero,
949
950 list.push_back(p1);
951 list.push_back(p2);
952 list.push_back(p3);
953 list.push_back(p4);
954 list.push_back(p5);
955 list.push_back(p6);
956 list.push_back(p7);
957 list.push_back(p8);
958 list.push_back(p9);
959 } else {
960 // Default condition
961 if(rdm < (1.-kaonicFSprob)){ // pion
962 INCL_DEBUG("pionic pp final stat
963 sum = read_file(dataPathppbar, p
964 rdm = (rdm/(1.-kaonicFSprob))*su
965 //now get the line number in the
966 G4int n = findStringNumber(rdm,
967 for(G4int j = 0; j < static_cast
968 if(particle_types[n][j] == "pi
969 Particle *p = new Particle
970 list.push_back(p);
971 }
972 else if(particle_types[n][j] =
973 Particle *p = new Particle
974 list.push_back(p);
975 }
976 else if(particle_types[n][j] =
977 Particle *p = new Particle
978 list.push_back(p);
979 }
980 else if(particle_types[n][j] =
981 Particle *p = new Particle
982 list.push_back(p);
983 }
984 else if(particle_types[n][j] =
985 Particle *p = new Particle
986 list.push_back(p);
987 }
988 else{
989 INCL_ERROR("Some non-existing FS
990 for(G4int jj = 0; jj < static_ca
991 std::cout << "gotcha! " << par
992 }
993 }
994 } //end of pionic option
995 else{
996 INCL_DEBUG("kaonic pp final stat
997 sum = read_file(dataPathppbark,
998 rdm = ((1-rdm)/kaonicFSprob)*sum
999 //now get the line number in the
1000 G4int n = findStringNumber(rdm,
1001 for(G4int j = 0; j < static_cast<
1002 if(particle_types[n][j] == "p
1003 Particle *p = new Particl
1004 list.push_back(p);
1005 }
1006 else if(particle_types[n][j]
1007 Particle *p = new Particl
1008 list.push_back(p);
1009 }
1010 else if(particle_types[n][j]
1011 Particle *p = new Particl
1012 list.push_back(p);
1013 }
1014 else if(particle_types[n][j]
1015 Particle *p = new Particl
1016 list.push_back(p);
1017 }
1018 else if(particle_types[n][j]
1019 Particle *p = new Particl
1020 list.push_back(p);
1021 }
1022 else if(particle_types[n][j]
1023 Particle *p = new Particl
1024 list.push_back(p);
1025 }
1026 else if(particle_types[n][j]
1027 Particle *p = new Particl
1028 list.push_back(p);
1029 }
1030 else if(particle_types[n][j]
1031 Particle *p = new Particl
1032 list.push_back(p);
1033 }
1034 else if(particle_types[n][j]
1035 Particle *p = new Particl
1036 list.push_back(p);
1037 }
1038 else{
1039 INCL_ERROR("Some non-exis
1040 for(G4int jj = 0; jj < st
1041 std::cout << "gotcha! "
1042 }
1043 }
1044 }
1045 } // end of kaonic option
1046 } // end of default condition
1047 } // end of ppbar and nnbar case
1048
1049
1050
1051 nucleon->setType(list[0]->getType());
1052 antinucleon->setType(list[1]->getType());
1053
1054 ParticleList finallist;
1055
1056 finallist.push_back(nucleon);
1057 finallist.push_back(antinucleon);
1058
1059 if(list.size() > 2){
1060 for (G4int i = 2; i < (G4int)(list.size
1061 finallist.push_back(list[i]);
1062 }
1063 }
1064
1065 if(finallist.size()==2){
1066 G4double mn=nucleon->getMass();
1067 G4double my=antinucleon->getMass();
1068
1069 G4double ey=(sqrtS*sqrtS+my*my-mn*mn)/(
1070 G4double en=std::sqrt(ey*ey-my*my+mn*mn
1071 nucleon->setEnergy(en);
1072 antinucleon->setEnergy(ey);
1073 G4double py=std::sqrt(ey*ey-my*my);
1074
1075 ThreeVector mom_antinucleon = Random::n
1076 antinucleon->setMomentum(mom_antinucleo
1077 nucleon->setMomentum(-mom_antinucleon);
1078 }
1079 else if(finallist.size() > 2){
1080 PhaseSpaceGenerator::generate(sqrtS, fi
1081 }
1082 else{
1083 INCL_ERROR("less than 2 mesons in NNbar
1084 }
1085
1086
1087 for (G4int i = 0; i < 2; i++) {
1088 fs->addModifiedParticle(finallist[i]);
1089 }
1090 if(finallist.size()>2){
1091 for (G4int i = 2; i < (G4int)(list.size
1092 fs->addCreatedParticle(finallist[i]);
1093 }
1094 }
1095
1096
1097 //fs->addDestroyedParticle(nucleon);
1098 //fs->addDestroyedParticle(antinucleon);
1099
1100
1101 }
1102 }
1103
1104