Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCLNNToMultiPionsChannel.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#include "G4INCLNNToMultiPionsChannel.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
47namespace G4INCL {
48
49 const G4double NNToMultiPionsChannel::angularSlope = 6.;
50
51 NNToMultiPionsChannel::NNToMultiPionsChannel(const G4int npi, Particle *p1, Particle *p2)
52 : npion(npi),
53 iso1(0),
54 iso2(0),
55 particle1(p1),
56 particle2(p2)
57 {
58 std::fill(isosp, isosp+4, 0);
59 }
60
64
65 void NNToMultiPionsChannel::fillFinalState(FinalState *fs) {
66// assert(npion > 0 && npion < 5);
67
68 iso1=ParticleTable::getIsospin(particle1->getType());
69 iso2=ParticleTable::getIsospin(particle2->getType());
70
71 ParticleList list;
72 list.push_back(particle1);
73 list.push_back(particle2);
74 fs->addModifiedParticle(particle1);
75 fs->addModifiedParticle(particle2);
76
77 isospinRepartition();
78
79 const ParticleType tn1=ParticleTable::getNucleonType(iso1);
80 particle1->setType(tn1);
81 const ParticleType tn2=ParticleTable::getNucleonType(iso2);
82 particle2->setType(tn2);
83 const ThreeVector &rcolnucleon1 = particle1->getPosition();
84 const ThreeVector &rcolnucleon2 = particle2->getPosition();
85 const ThreeVector rcol = (rcolnucleon1+rcolnucleon2)*0.5;
86 const ThreeVector zero;
87 for(G4int i=0; i<npion; ++i) {
88 const ParticleType pionType=ParticleTable::getPionType(isosp[i]);
89 Particle *pion = new Particle(pionType,zero,rcol);
90 list.push_back(pion);
91 fs->addCreatedParticle(pion);
92 }
93
94 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2);
95 G4int biasIndex = ((Random::shoot()<0.5) ? 0 : 1);
96 PhaseSpaceGenerator::generateBiased(sqrtS, list, biasIndex, angularSlope);
97
98 }
99
100 void NNToMultiPionsChannel::isospinRepartition() {
101 const G4double rjcd=Random::shoot();
102 G4double p;
103 const G4int itot=iso1+iso2;
104
105 if (npion == 1) {
106 p=3.*rjcd;
107 if (p < 1.) pn_ppPim();
108 else if (p < 2.) pn_pnPi0();
109 else pn_nnPip();
110 }
111 else if (npion == 2) {
112 if (itot == 2) {
113 p=20.*rjcd;
114 if (p >= 14.) pp_nnPipPip();
115 else if (p >= 11.) pp_pnPipPi0();
116 else if (p >= 7.) pp_ppPi0Pi0();
117 else pp_ppPipPim();
118 }
119 else if (itot == -2) {
120 p=20.*rjcd;
121 if (p >= 14.) nn_ppPimPim();
122 else if (p >= 11.) nn_pnPimPi0();
123 else if (p >= 7.) nn_nnPi0Pi0();
124 else nn_nnPipPim();
125 }
126 else {
127 G4double pp=Random::shoot();
128 if (pp > 0.5) {
129 p=3.*rjcd;
130 if (p < 2.) {
131 pn_pnPipPim();
132 }
133 else {
134 pn_pnPi0Pi0();
135 }
136 }
137 else {
138 p=60.*rjcd;
139 if (p >= 51.) pn_nnPipPi0();
140 else if (p >= 33.) pn_pnPi0Pi0();
141 else if (p >= 9.) pn_pnPipPim();
142 else pn_ppPimPi0();
143 }
144 }
145 }
146 else if (npion == 3) {
147 p=60.*rjcd;
148 if (itot == 2) {
149 if (p >= 42.) pp_nnPipPipPi0();
150 else if (p >= 39.) pp_pnPipPi0Pi0();
151 else if (p >= 33.) pp_pnPipPipPim();
152 else if (p >= 22.) pp_ppPi0Pi0Pi0();
153 else pp_ppPipPimPi0();
154 }
155 else if (itot == -2) {
156 if (p >= 42.) nn_ppPimPimPi0();
157 else if (p >= 39.) nn_pnPimPi0Pi0();
158 else if (p >= 33.) nn_pnPipPimPim();
159 else if (p >= 22.) nn_nnPi0Pi0Pi0();
160 else nn_nnPipPimPi0();
161 }
162 else {
163 if (p >= 57.) pn_nnPipPi0Pi0();
164 else if (p >= 51.) pn_nnPipPipPim();
165 else if (p >= 37.) pn_pnPi0Pi0Pi0();
166 else if (p >= 9.) pn_pnPi0PipPim();
167 else if (p >= 6.) pn_ppPimPi0Pi0();
168 else pn_ppPimPimPip();
169
170 }
171 }
172 else if (npion == 4) {
173 p=60.*rjcd;
174 if (itot == 2) {
175 if (p >= 48.) pp_nnPipPipPipPim();
176 else if (p >= 42.) pp_nnPipPipPi0Pi0();
177 else if (p >= 36.) pp_pnPipPipPi0Pim();
178 else if (p >= 33.) pp_pnPipPi0Pi0Pi0();
179 else if (p >= 19.) pp_ppPipPipPimPim();
180 else if (p >= 4.) pp_ppPipPi0Pi0Pim();
181 else pp_ppPi0Pi0Pi0Pi0();
182 }
183 else if (itot == -2) {
184 if (p >= 48.) nn_ppPipPimPimPim();
185 else if (p >= 42.) nn_ppPi0Pi0PimPim();
186 else if (p >= 36.) nn_pnPipPi0PimPim();
187 else if (p >= 33.) nn_pnPi0Pi0Pi0Pim();
188 else if (p >= 19.) nn_nnPipPipPimPim();
189 else if (p >= 4.) nn_nnPipPi0Pi0Pim();
190 else nn_nnPi0Pi0Pi0Pi0();
191 }
192 else {
193 G4double pp=Random::shoot();
194 if (pp > 0.5) {
195 p=9.*rjcd;
196 if (p < 1.) pn_pnPi0Pi0Pi0Pi0();
197 else if (p < 5.) pn_pnPipPi0Pi0Pim();
198 else pn_pnPipPipPimPim();
199 }
200 else {
201 if (p < 3.) pn_ppPi0Pi0Pi0Pim();
202 else if (p < 9.) pn_ppPipPi0PimPim();
203 else if (p < 15.) pn_pnPi0Pi0Pi0Pi0();
204 else if (p < 35.) pn_pnPipPi0Pi0Pim();
205 else if (p < 51.) pn_pnPipPipPimPim();
206 else if (p < 54.) pn_nnPipPi0Pi0Pi0();
207 else pn_nnPipPipPi0Pim();
208 }
209 }
210 }
211
212 std::shuffle(isosp,isosp+npion,Random::getAdapter());
213 inter2Part(0.5);
214 }
215
216
217 void NNToMultiPionsChannel::pn_ppPim() {
218 isosp[0]=-2;
219 iso1=1;
220 iso2=1;
221 }
222 void NNToMultiPionsChannel::pn_pnPi0() {
223 isosp[0]=0;
224 }
225 void NNToMultiPionsChannel::pn_nnPip() {
226 isosp[0]=2;
227 iso1=-1;
228 iso2=-1;
229 }
230 void NNToMultiPionsChannel::pp_nnPipPip() {
231 isosp[0]=2;
232 isosp[1]=2;
233 iso1=-1;
234 iso2=-1;
235 }
236 void NNToMultiPionsChannel::nn_ppPimPim() {
237 isosp[0]=-2;
238 isosp[1]=-2;
239 iso1=1;
240 iso2=1;
241 }
242 void NNToMultiPionsChannel::pn_pnPipPim() {
243 isosp[0]=2;
244 isosp[1]=-2;
245 }
246 void NNToMultiPionsChannel::pn_pnPi0Pi0() {
247 isosp[0]=0;
248 isosp[1]=0;
249 }
250 void NNToMultiPionsChannel::pp_ppPipPim() {
251 isosp[0]=2;
252 isosp[1]=-2;
253 }
254 void NNToMultiPionsChannel::nn_nnPipPim() {
255 isosp[0]=2;
256 isosp[1]=-2;
257 }
258 void NNToMultiPionsChannel::pp_ppPi0Pi0() {
259 isosp[0]=0;
260 isosp[1]=0;
261 }
262 void NNToMultiPionsChannel::nn_nnPi0Pi0() {
263 isosp[0]=0;
264 isosp[1]=0;
265 }
266 void NNToMultiPionsChannel::pp_pnPipPi0() {
267 isosp[0]=2;
268 isosp[1]=0;
269 iso1=1;
270 iso2=-1;
271 }
272 void NNToMultiPionsChannel::pn_ppPimPi0() {
273 isosp[0]=-2;
274 isosp[1]=0;
275 iso1=1;
276 iso2=1;
277 }
278 void NNToMultiPionsChannel::pn_nnPipPi0() {
279 isosp[0]=2;
280 isosp[1]=0;
281 iso1=-1;
282 iso2=-1;
283 }
284 void NNToMultiPionsChannel::nn_pnPimPi0() {
285 isosp[0]=-2;
286 isosp[1]=0;
287 iso1=1;
288 iso2=-1;
289 }
290 void NNToMultiPionsChannel::pp_pnPipPi0Pi0() {
291 isosp[0]=2;
292 isosp[1]=0;
293 isosp[2]=0;
294 iso1=1;
295 iso2=-1;
296 }
297 void NNToMultiPionsChannel::nn_pnPimPi0Pi0() {
298 isosp[0]=-2;
299 isosp[1]=0;
300 isosp[2]=0;
301 iso1=1;
302 iso2=-1;
303 }
304 void NNToMultiPionsChannel::pn_nnPipPi0Pi0() {
305 isosp[0]=2;
306 isosp[1]=0;
307 isosp[2]=0;
308 iso1=-1;
309 iso2=-1;
310 }
311 void NNToMultiPionsChannel::pp_ppPipPimPi0() {
312 isosp[0]=2;
313 isosp[1]=-2;
314 isosp[2]=0;
315 }
316 void NNToMultiPionsChannel::nn_nnPipPimPi0() {
317 isosp[0]=2;
318 isosp[1]=-2;
319 isosp[2]=0;
320 }
321 void NNToMultiPionsChannel::pp_ppPi0Pi0Pi0() {
322 isosp[0]=0;
323 isosp[1]=0;
324 isosp[2]=0;
325 }
326 void NNToMultiPionsChannel::nn_nnPi0Pi0Pi0() {
327 isosp[0]=0;
328 isosp[1]=0;
329 isosp[2]=0;
330 }
331 void NNToMultiPionsChannel::pp_pnPipPipPim() {
332 isosp[0]=2;
333 isosp[1]=2;
334 isosp[2]=-2;
335 iso1=1;
336 iso2=-1;
337 }
338 void NNToMultiPionsChannel::pp_nnPipPipPi0() {
339 isosp[0]=2;
340 isosp[1]=2;
341 isosp[2]=0;
342 iso1=-1;
343 iso2=-1;
344 }
345 void NNToMultiPionsChannel::pn_ppPimPi0Pi0() {
346 isosp[0]=-2;
347 isosp[1]=0;
348 isosp[2]=0;
349 iso1=1;
350 iso2=1;
351 }
352 void NNToMultiPionsChannel::pn_ppPimPimPip() {
353 isosp[0]=-2;
354 isosp[1]=-2;
355 isosp[2]=2;
356 iso1=1;
357 iso2=1;
358 }
359 void NNToMultiPionsChannel::pn_pnPi0PipPim() {
360 isosp[0]=0;
361 isosp[1]=2;
362 isosp[2]=-2;
363 }
364 void NNToMultiPionsChannel::pn_pnPi0Pi0Pi0() {
365 isosp[0]=0;
366 isosp[1]=0;
367 isosp[2]=0;
368 }
369 void NNToMultiPionsChannel::pn_nnPipPipPim() {
370 isosp[0]=2;
371 isosp[1]=2;
372 isosp[2]=-2;
373 iso1=-1;
374 iso2=-1;
375 }
376 void NNToMultiPionsChannel::nn_pnPipPimPim() {
377 isosp[0]=2;
378 isosp[1]=-2;
379 isosp[2]=-2;
380 iso1=1;
381 iso2=-1;
382 }
383 void NNToMultiPionsChannel::nn_ppPimPimPi0() {
384 isosp[0]=-2;
385 isosp[1]=-2;
386 isosp[2]=0;
387 iso1=1;
388 iso2=1;
389 }
390 void NNToMultiPionsChannel::pp_nnPipPipPi0Pi0() {
391 isosp[0]=2;
392 isosp[1]=2;
393 isosp[2]=0;
394 isosp[3]=0;
395 iso1=-1;
396 iso2=-1;
397 }
398 void NNToMultiPionsChannel::pp_nnPipPipPipPim() {
399 isosp[0]=2;
400 isosp[1]=2;
401 isosp[2]=2;
402 isosp[3]=-2;
403 iso1=-1;
404 iso2=-1;
405 }
406 void NNToMultiPionsChannel::nn_ppPi0Pi0PimPim() {
407 isosp[0]=0;
408 isosp[1]=0;
409 isosp[2]=-2;
410 isosp[3]=-2;
411 iso1=1;
412 iso2=1;
413 }
414 void NNToMultiPionsChannel::nn_ppPipPimPimPim() {
415 isosp[0]=2;
416 isosp[1]=-2;
417 isosp[2]=-2;
418 isosp[3]=-2;
419 iso1=1;
420 iso2=1;
421 }
422 void NNToMultiPionsChannel::pp_ppPi0Pi0Pi0Pi0() {
423 isosp[0]=0;
424 isosp[1]=0;
425 isosp[2]=0;
426 isosp[3]=0;
427 }
428 void NNToMultiPionsChannel::nn_nnPi0Pi0Pi0Pi0() {
429 isosp[0]=0;
430 isosp[1]=0;
431 isosp[2]=0;
432 isosp[3]=0;
433 }
434 void NNToMultiPionsChannel::pn_pnPi0Pi0Pi0Pi0() {
435 isosp[0]=0;
436 isosp[1]=0;
437 isosp[2]=0;
438 isosp[3]=0;
439 }
440 void NNToMultiPionsChannel::pp_ppPipPi0Pi0Pim() {
441 isosp[0]=2;
442 isosp[1]=0;
443 isosp[2]=0;
444 isosp[3]=-2;
445 }
446 void NNToMultiPionsChannel::nn_nnPipPi0Pi0Pim() {
447 isosp[0]=2;
448 isosp[1]=0;
449 isosp[2]=0;
450 isosp[3]=-2;
451 }
452 void NNToMultiPionsChannel::pn_pnPipPi0Pi0Pim() {
453 isosp[0]=2;
454 isosp[1]=0;
455 isosp[2]=0;
456 isosp[3]=-2;
457 }
458 void NNToMultiPionsChannel::pp_ppPipPipPimPim() {
459 isosp[0]=2;
460 isosp[1]=2;
461 isosp[2]=-2;
462 isosp[3]=-2;
463 }
464 void NNToMultiPionsChannel::nn_nnPipPipPimPim() {
465 isosp[0]=2;
466 isosp[1]=2;
467 isosp[2]=-2;
468 isosp[3]=-2;
469 }
470 void NNToMultiPionsChannel::pn_pnPipPipPimPim() {
471 isosp[0]=2;
472 isosp[1]=2;
473 isosp[2]=-2;
474 isosp[3]=-2;
475 }
476 void NNToMultiPionsChannel::pp_pnPipPi0Pi0Pi0() {
477 isosp[0]=2;
478 isosp[1]=0;
479 isosp[2]=0;
480 isosp[3]=0;
481 iso1=1;
482 iso2=-1;
483 }
484 void NNToMultiPionsChannel::pn_nnPipPi0Pi0Pi0() {
485 isosp[0]=2;
486 isosp[1]=0;
487 isosp[2]=0;
488 isosp[3]=0;
489 iso1=-1;
490 iso2=-1;
491 }
492 void NNToMultiPionsChannel::pp_nnPipPi0Pi0Pi0() {
493 isosp[0]=2;
494 isosp[1]=0;
495 isosp[2]=0;
496 isosp[3]=0;
497 iso1=-1;
498 iso2=-1;
499 }
500 void NNToMultiPionsChannel::pp_pnPipPipPi0Pim() {
501 isosp[0]=2;
502 isosp[1]=2;
503 isosp[2]=0;
504 isosp[3]=-2;
505 iso1=1;
506 iso2=-1;
507 }
508 void NNToMultiPionsChannel::pn_nnPipPipPi0Pim() {
509 isosp[0]=2;
510 isosp[1]=2;
511 isosp[2]=0;
512 isosp[3]=-2;
513 iso1=-1;
514 iso2=-1;
515 }
516 void NNToMultiPionsChannel::pp_nnPipPipPi0Pim() {
517 isosp[0]=2;
518 isosp[1]=2;
519 isosp[2]=0;
520 isosp[3]=-2;
521 iso1=-1;
522 iso2=-1;
523 }
524 void NNToMultiPionsChannel::nn_pnPi0Pi0Pi0Pim() {
525 isosp[0]=0;
526 isosp[1]=0;
527 isosp[2]=0;
528 isosp[3]=-2;
529 iso1=1;
530 iso2=-1;
531 }
532 void NNToMultiPionsChannel::pn_ppPi0Pi0Pi0Pim() {
533 isosp[0]=0;
534 isosp[1]=0;
535 isosp[2]=0;
536 isosp[3]=-2;
537 iso1=1;
538 iso2=1;
539 }
540 void NNToMultiPionsChannel::nn_pnPipPi0PimPim() {
541 isosp[0]=2;
542 isosp[1]=0;
543 isosp[2]=-2;
544 isosp[3]=-2;
545 iso1=1;
546 iso2=-1;
547 }
548 void NNToMultiPionsChannel::pn_ppPipPi0PimPim() {
549 isosp[0]=2;
550 isosp[1]=0;
551 isosp[2]=-2;
552 isosp[3]=-2;
553 iso1=1;
554 iso2=1;
555 }
556
557 void NNToMultiPionsChannel::inter2Part(const G4double p) {
558
559 if (Random::shoot() < p) std::swap(iso1,iso2);
560
561 }
562
563
564}
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4INCL::FinalState::addModifiedParticle
void addModifiedParticle(Particle *p)
Definition G4INCLFinalState.cc:60
G4INCL::FinalState::addCreatedParticle
void addCreatedParticle(Particle *p)
Definition G4INCLFinalState.cc:75
G4INCL::NNToMultiPionsChannel::NNToMultiPionsChannel
NNToMultiPionsChannel(const G4int, Particle *, Particle *)
Definition G4INCLNNToMultiPionsChannel.cc:51
G4INCL::Particle::getPosition
const G4INCL::ThreeVector & getPosition() const
Definition G4INCLParticle.hh:956
G4INCL::Particle::getType
G4INCL::ParticleType getType() const
Definition G4INCLParticle.hh:187
G4INCL::Particle::setType
void setType(ParticleType t)
Definition G4INCLParticle.hh:196
G4INCL::KinematicsUtils::totalEnergyInCM
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
Definition G4INCLKinematicsUtils.cc:121
G4INCL::ParticleTable::getIsospin
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
Definition G4INCLParticleTable.cc:531
G4INCL::ParticleTable::getNucleonType
ParticleType getNucleonType(const G4int isosp)
Get the type of nucleon.
Definition G4INCLParticleTable.cc:1455
G4INCL::ParticleTable::getPionType
ParticleType getPionType(const G4int isosp)
Get the type of pion.
Definition G4INCLParticleTable.cc:1442
G4INCL::PhaseSpaceGenerator::generateBiased
void generateBiased(const G4double sqrtS, ParticleList &particles, const size_t index, const G4double slope)
Generate a biased event in the CM system.
Definition G4INCLPhaseSpaceGenerator.cc:98
G4INCL::Random::getAdapter
Adapter const & getAdapter()
Definition G4INCLRandom.cc:257
G4INCL::Random::shoot
G4double shoot()
Definition G4INCLRandom.cc:93