| 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
38 /*
39 * G4INCLBinaryCollisionAvatar.cc
40 *
41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi
43 */
44
45 #include "G4INCLBinaryCollisionAvatar.hh"
46 #include "G4INCLElasticChannel.hh"
47 #include "G4INCLRecombinationChannel.hh"
48 #include "G4INCLDeltaProductionChannel.hh"
49 #include "G4INCLNNToMultiPionsChannel.hh"
50 #include "G4INCLNNToNNEtaChannel.hh"
51 #include "G4INCLNDeltaEtaProductionChannel.hh"
52 #include "G4INCLNNEtaToMultiPionsChannel.hh"
53 #include "G4INCLNNToNNOmegaChannel.hh"
54 #include "G4INCLNDeltaOmegaProductionChannel.h
55 #include "G4INCLNNOmegaToMultiPionsChannel.hh"
56 #include "G4INCLCrossSections.hh"
57 #include "G4INCLKinematicsUtils.hh"
58 #include "G4INCLRandom.hh"
59 #include "G4INCLParticleTable.hh"
60 #include "G4INCLPauliBlocking.hh"
61 #include "G4INCLPiNElasticChannel.hh"
62 #include "G4INCLPiNToDeltaChannel.hh"
63 #include "G4INCLPiNToMultiPionsChannel.hh"
64 #include "G4INCLPiNToEtaChannel.hh"
65 #include "G4INCLPiNToOmegaChannel.hh"
66 #include "G4INCLEtaNElasticChannel.hh"
67 #include "G4INCLEtaNToPiNChannel.hh"
68 #include "G4INCLEtaNToPiPiNChannel.hh"
69 #include "G4INCLOmegaNElasticChannel.hh"
70 #include "G4INCLOmegaNToPiNChannel.hh"
71 #include "G4INCLNNToNLKChannel.hh"
72 #include "G4INCLNNToNSKChannel.hh"
73 #include "G4INCLNNToNLKpiChannel.hh"
74 #include "G4INCLNNToNSKpiChannel.hh"
75 #include "G4INCLNNToNLK2piChannel.hh"
76 #include "G4INCLNNToNSK2piChannel.hh"
77 #include "G4INCLNNToNNKKbChannel.hh"
78 #include "G4INCLNNToMissingStrangenessChannel.
79 #include "G4INCLNDeltaToNLKChannel.hh"
80 #include "G4INCLNDeltaToNSKChannel.hh"
81 #include "G4INCLNDeltaToDeltaLKChannel.hh"
82 #include "G4INCLNDeltaToDeltaSKChannel.hh"
83 #include "G4INCLNDeltaToNNKKbChannel.hh"
84 #include "G4INCLNpiToLKChannel.hh"
85 #include "G4INCLNpiToSKChannel.hh"
86 #include "G4INCLNpiToLKpiChannel.hh"
87 #include "G4INCLNpiToSKpiChannel.hh"
88 #include "G4INCLNpiToLK2piChannel.hh"
89 #include "G4INCLNpiToSK2piChannel.hh"
90 #include "G4INCLNpiToNKKbChannel.hh"
91 #include "G4INCLNpiToMissingStrangenessChannel
92 #include "G4INCLNKElasticChannel.hh"
93 #include "G4INCLNKToNKChannel.hh"
94 #include "G4INCLNKToNKpiChannel.hh"
95 #include "G4INCLNKToNK2piChannel.hh"
96 #include "G4INCLNKbElasticChannel.hh"
97 #include "G4INCLNKbToNKbChannel.hh"
98 #include "G4INCLNKbToNKbpiChannel.hh"
99 #include "G4INCLNKbToNKb2piChannel.hh"
100 #include "G4INCLNKbToLpiChannel.hh"
101 #include "G4INCLNKbToL2piChannel.hh"
102 #include "G4INCLNKbToSpiChannel.hh"
103 #include "G4INCLNKbToS2piChannel.hh"
104 #include "G4INCLNYElasticChannel.hh"
105 #include "G4INCLNLToNSChannel.hh"
106 #include "G4INCLNSToNLChannel.hh"
107 #include "G4INCLNSToNSChannel.hh"
108 #include "G4INCLOmegaNToPiPiNChannel.hh"
109 #include "G4INCLStore.hh"
110 #include "G4INCLBook.hh"
111 #include "G4INCLLogger.hh"
112 #include <string>
113 #include <sstream>
114 // #include <cassert>
115 #include "G4INCLNNbarElasticChannel.hh"
116 #include "G4INCLNNbarCEXChannel.hh"
117 #include "G4INCLNNbarToLLbarChannel.hh"
118 #include "G4INCLNNbarToNNbarpiChannel.hh"
119 #include "G4INCLNNbarToNNbar2piChannel.hh"
120 #include "G4INCLNNbarToNNbar3piChannel.hh"
121 #include "G4INCLNNbarToAnnihilationChannel.hh"
122
123 namespace G4INCL {
124
125 // WARNING: if you update the default cutNN
126 // cutNNSquared variable, too.
127 G4ThreadLocal G4double BinaryCollisionAvatar
128 G4ThreadLocal G4double BinaryCollisionAvatar
129 G4ThreadLocal G4double BinaryCollisionAvatar
130
131 BinaryCollisionAvatar::BinaryCollisionAvatar
132 G4INCL::Nucleus *n, G4INCL::Particle *p1
133 : InteractionAvatar(time, n, p1, p2), theC
134 isParticle1Spectator(false),
135 isParticle2Spectator(false),
136 isElastic(false),
137 isStrangeProduction(false)
138 {
139 setType(CollisionAvatarType);
140 }
141
142 BinaryCollisionAvatar::~BinaryCollisionAvata
143 }
144
145 G4INCL::IChannel* BinaryCollisionAvatar::get
146 // We already check cutNN at avatar creati
147 // again here. For composite projectiles,
148 // avatars with no cutNN before any collis
149 if(particle1->isNucleon()
150 && particle2->isNucleon()
151 && theNucleus->getStore()->getBook().g
152 const G4double energyCM2 = KinematicsUti
153 // Below a certain cut value we don't do
154 if(energyCM2 < cutNNSquared) {
155 INCL_DEBUG("CM energy = sqrt(" << ener
156 << ") MeV = cutNN" << "; returning
157 InteractionAvatar::restoreParticles();
158 return NULL;
159 }
160 }
161
162 /** Check again the distance of approach.
163 * realised, we have to perform a check in
164 * distance four-vector as
165 * \f[ (0, \Delta\vec{x}), \f]
166 * where \f$\Delta\vec{x}\f$ is the distan
167 * their minimum distance of approach (i.e
168 * boosting this four-vector to the CM fra
169 * obtain a new four vector
170 * \f[ (\Delta t', \Delta\vec{x}'), \f]
171 * with a non-zero time component (the col
172 * the two particles in the lab system, bu
173 * for the avatar to be realised, we requi
174 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/
175 * Note that \f$|\Delta\vec{x}'|\leq|\Delt
176 * above is more restrictive than the chec
177 * G4INCL::Propagation::StandardPropagatio
178 * In other words, the avatar generation c
179 * avatars.
180 */
181 ThreeVector minimumDistance = particle1->g
182 minimumDistance -= particle2->getPosition(
183 const G4double betaDotX = boostVector.dot(
184 const G4double minDist = Math::tenPi*(mini
185 if(minDist > theCrossSection) {
186 INCL_DEBUG("CM distance of approach is t
187 theCrossSection <<"; returning a NULL
188 InteractionAvatar::restoreParticles();
189 return NULL;
190 }
191
192 /** Bias apply for this reaction in order
193 * ParticleBias for all stange particles.
194 * Can be reduced after because of the saf
195 */
196 G4double bias_apply = 1.;
197 if(bias != 1.) bias_apply = Particle::getB
198
199 //// NN
200 if(particle1->isNucleon() && particle2->is
201
202 G4double NLKProductionCX = CrossSections
203 G4double NSKProductionCX = CrossSections
204 G4double NLKpiProductionCX = CrossSectio
205 G4double NSKpiProductionCX = CrossSectio
206 G4double NLK2piProductionCX = CrossSecti
207 G4double NSK2piProductionCX = CrossSecti
208 G4double NNKKbProductionCX = CrossSectio
209 G4double NNMissingCX = CrossSections::NN
210
211 const G4double UnStrangeProdCX = CrossSe
212 + CrossSect
213 + CrossSect
214 + CrossSect
215 + CrossSect
216 + CrossSect
217 const G4double StrangenessProdCX = (NLKP
218
219 G4double counterweight = (1. - bias_appl
220
221 if(counterweight < 0.5) {
222 counterweight = 0.5;
223 bias_apply = 0.5*UnStrangeProdCX/Stra
224 NLKProductionCX = CrossSections::NNTo
225 NSKProductionCX = CrossSections::NNTo
226 NLKpiProductionCX = CrossSections::NN
227 NSKpiProductionCX = CrossSections::NN
228 NLK2piProductionCX = CrossSections::N
229 NSK2piProductionCX = CrossSections::N
230 NNKKbProductionCX = CrossSections::NN
231 NNMissingCX = CrossSections::NNToMiss
232 }
233
234
235 const G4double elasticCX = CrossSections
236 const G4double deltaProductionCX = Cross
237 const G4double onePiProductionCX = Cross
238 const G4double twoPiProductionCX = Cross
239 const G4double threePiProductionCX = Cro
240 const G4double fourPiProductionCX = Cros
241
242 const G4double etaProductionCX = CrossSe
243 const G4double etadeltaProductionCX = Cr
244 const G4double etaonePiProductionCX = Cr
245 const G4double etatwoPiProductionCX = Cr
246 const G4double etathreePiProductionCX =
247 const G4double etafourPiProductionCX = C
248 const G4double omegaProductionCX = Cross
249 const G4double omegadeltaProductionCX =
250 const G4double omegaonePiProductionCX =
251 const G4double omegatwoPiProductionCX =
252 const G4double omegathreePiProductionCX
253 const G4double omegafourPiProductionCX =
254
255 const G4double totCX=CrossSections::tota
256
257 // assert(std::fabs(totCX-elasticCX-deltaProdu
258
259 const G4double rChannel=Random::shoot()
260
261 if(elasticCX > rChannel) {
262 // Elastic NN channel
263 isElastic = true;
264 INCL_DEBUG("NN interaction: elastic ch
265 weight = counterweight;
266 return new ElasticChannel(particle1, p
267 } else if((elasticCX + deltaProductionCX
268 isElastic = false;
269 // NN -> N Delta channel is chosen
270 INCL_DEBUG("NN interaction: Delta chan
271 weight = counterweight;
272 return new DeltaProductionChannel(part
273 } else if(elasticCX + deltaProductionCX
274 isElastic = false;
275 // NN -> PiNN channel is chosen
276 INCL_DEBUG("NN interaction: one Pion c
277 weight = counterweight;
278 return new NNToMultiPionsChannel(1,par
279 } else if(elasticCX + deltaProductionCX
280 isElastic = false;
281 // NN -> 2PiNN channel is chosen
282 INCL_DEBUG("NN interaction: two Pions
283 weight = counterweight;
284 return new NNToMultiPionsChannel(2,par
285 } else if(elasticCX + deltaProductionCX
286 isElastic = false;
287 // NN -> 3PiNN channel is chosen
288 INCL_DEBUG("NN interaction: three Pion
289 weight = counterweight;
290 return new NNToMultiPionsChannel(3,par
291 } else if(elasticCX + deltaProductionCX
292 isElastic = false;
293 // NN -> 4PiNN channel is chosen
294 INCL_DEBUG("NN interaction: four
295 weight = counterweight;
296 return new NNToMultiPionsChannel
297 } else if(elasticCX + deltaProductionCX
298 + etaProductionCX >
299 isElastic = false;
300 // NN -> NNEta channel is chosen
301 INCL_DEBUG("NN interaction: Eta channel
302 weight = counterweight;
303 return new NNToNNEtaChannel(particle1,
304 } else if(elasticCX + deltaProductionCX
305 + etaProductionCX +
306 isElastic = false;
307 // NN -> N Delta Eta channel is chosen
308 INCL_DEBUG("NN interaction: Delta Eta c
309 weight = counterweight;
310 return new NDeltaEtaProductionChannel(p
311 } else if(elasticCX + deltaProductionCX
312 + etaProductionCX +
313 isElastic = false;
314 // NN -> EtaPiNN channel is chosen
315 INCL_DEBUG("NN interaction: Eta + one P
316 weight = counterweight;
317 return new NNEtaToMultiPionsChannel(1,p
318 } else if(elasticCX + deltaProductionCX
319 + etaProductionCX +
320 isElastic = false;
321 // NN -> Eta2PiNN channel is chosen
322 INCL_DEBUG("NN interaction: Eta + two P
323 weight = counterweight;
324 return new NNEtaToMultiPionsChannel(2,p
325 } else if(elasticCX + deltaProductionCX
326 + etaProductionCX +
327 isElastic = false;
328 // NN -> Eta3PiNN channel is chosen
329 INCL_DEBUG("NN interaction: Eta + three
330 weight = counterweight;
331 return new NNEtaToMultiPionsChannel(3,p
332 } else if(elasticCX + deltaProductionCX
333 + etaProductionCX +
334 isElastic = false;
335 // NN -> Eta4PiNN channel is chosen
336 INCL_DEBUG("NN interaction: Eta + four
337 weight = counterweight;
338 return new NNEtaToMultiPionsChannel(4,p
339 } else if(elasticCX + deltaProductionCX
340 + etaProductionCX +
341 + omegaProductionCX
342 isElastic = false;
343 // NN -> NNOmega channel is chosen
344 INCL_DEBUG("NN interaction: Omega chann
345 weight = counterweight;
346 return new NNToNNOmegaChannel(particle1
347 } else if(elasticCX + deltaProductionCX
348 + etaProductionCX +
349 + omegaProductionCX
350 isElastic = false;
351 // NN -> N Delta Omega channel is chosen
352 INCL_DEBUG("NN interaction: Delta Omega
353 weight = counterweight;
354 return new NDeltaOmegaProductionChannel
355 } else if(elasticCX + deltaProductionCX
356 + etaProductionCX +
357 + omegaProductionCX
358 isElastic = false;
359 // NN -> OmegaPiNN channel is chosen
360 INCL_DEBUG("NN interaction: Omega + one
361 weight = counterweight;
362 return new NNOmegaToMultiPionsChannel(1
363 } else if(elasticCX + deltaProductionCX
364 + etaProductionCX +
365 + omegaProductionCX
366 isElastic = false;
367 // NN -> Omega2PiNN channel is chosen
368 INCL_DEBUG("NN interaction: Omega + two
369 weight = counterweight;
370 return new NNOmegaToMultiPionsChannel(2
371 } else if(elasticCX + deltaProductionCX
372 + etaProductionCX +
373 + omegaProductionCX
374 isElastic = false;
375 // NN -> Omega3PiNN channel is chosen
376 INCL_DEBUG("NN interaction: Omega + thr
377 weight = counterweight;
378 return new NNOmegaToMultiPionsChannel(3
379 } else if(elasticCX + deltaProductionCX
380 + etaProductionCX +
381 + omegaProductionCX
382 isElastic = false;
383 // NN -> Omega4PiNN channel is chosen
384 INCL_DEBUG("NN interaction: Omega + fou
385 weight = counterweight;
386 return new NNOmegaToMultiPionsChannel(4
387 } else if(elasticCX + deltaProductionCX
388 + etaProductionCX +
389 + omegaProductionCX
390 + NLKProductionCX >
391 isElastic = false;
392 isStrangeProduction = true;
393 // NN -> NLK channel is chosen
394 INCL_DEBUG("NN interaction: NLK channe
395 weight = bias_apply;
396 return new NNToNLKChannel(particle1, p
397 } else if(elasticCX + deltaProductionCX
398 + etaProductionCX +
399 + omegaProductionCX
400 + NLKProductionCX +
401 isElastic = false;
402 isStrangeProduction = true;
403 // NN -> NLKpi channel is chosen
404 INCL_DEBUG("NN interaction: NLKpi chan
405 weight = bias_apply;
406 return new NNToNLKpiChannel(particle1,
407 } else if(elasticCX + deltaProductionCX
408 + etaProductionCX +
409 + omegaProductionCX
410 + NLKProductionCX +
411 isElastic = false;
412 isStrangeProduction = true;
413 // NN -> NLK2pi channel is chosen
414 INCL_DEBUG("NN interaction: NLK2pi cha
415 weight = bias_apply;
416 return new NNToNLK2piChannel(particle1
417 } else if(elasticCX + deltaProductionCX
418 + etaProductionCX +
419 + omegaProductionCX
420 + NLKProductionCX +
421 isElastic = false;
422 isStrangeProduction = true;
423 // NN -> NSK channel is chosen
424 INCL_DEBUG("NN interaction: NSK channe
425 weight = bias_apply;
426 return new NNToNSKChannel(particle1, p
427 } else if(elasticCX + deltaProductionCX
428 + etaProductionCX +
429 + omegaProductionCX
430 + NLKProductionCX +
431 isElastic = false;
432 isStrangeProduction = true;
433 // NN -> NSKpi channel is chosen
434 INCL_DEBUG("NN interaction: NSKpi chan
435 weight = bias_apply;
436 return new NNToNSKpiChannel(particle1,
437 } else if(elasticCX + deltaProductionCX
438 + etaProductionCX +
439 + omegaProductionCX
440 + NLKProductionCX +
441 isElastic = false;
442 isStrangeProduction = true;
443 // NN -> NSK2pi channel is chosen
444 INCL_DEBUG("NN interaction: NSK2pi cha
445 weight = bias_apply;
446 return new NNToNSK2piChannel(particle1
447 } else if(elasticCX + deltaProductionCX
448 + etaProductionCX +
449 + omegaProductionCX
450 + NLKProductionCX +
451 isElastic = false;
452 isStrangeProduction = true;
453 // NN -> NNKKb channel is chosen
454 INCL_DEBUG("NN interaction: NNKKb chan
455 weight = bias_apply;
456 return new NNToNNKKbChannel(particle1,
457 } else if(elasticCX + deltaProductionCX
458 + etaProductionCX +
459 + omegaProductionCX
460 + NLKProductionCX +
461 isElastic = false;
462 isStrangeProduction = true;
463 // NN -> Missing Strangeness channel is chosen
464 INCL_DEBUG("NN interaction: Missing St
465 weight = bias_apply;
466 return new NNToMissingStrangenessChann
467 } else {
468 INCL_WARN("inconsistency within the NN
469 if(NNMissingCX>0.) {
470 INCL_WARN("Returning an Missing St
471 weight = bias_apply;
472 isElastic = false;
473 isStrangeProduction = true;
474 return new NNToNNKKbChannel(partic
475 } else if(NNKKbProductionCX>0.) {
476 INCL_WARN("Returning an NNKKb chan
477 weight = bias_apply;
478 isElastic = false;
479 isStrangeProduction = true;
480 return new NNToNNKKbChannel(partic
481 } else if(NSK2piProductionCX>0.) {
482 INCL_WARN("Returning an NSK2pi cha
483 weight = bias_apply;
484 isElastic = false;
485 isStrangeProduction = true;
486 return new NNToNSK2piChannel(parti
487 } else if(NSKpiProductionCX>0.) {
488 INCL_WARN("Returning an NSKpi chan
489 weight = bias_apply;
490 isElastic = false;
491 isStrangeProduction = true;
492 return new NNToNSKpiChannel(partic
493 } else if(NSKProductionCX>0.) {
494 INCL_WARN("Returning an NSK channe
495 weight = bias_apply;
496 isElastic = false;
497 isStrangeProduction = true;
498 return new NNToNSKChannel(particle
499 } else if(NLK2piProductionCX>0.) {
500 INCL_WARN("Returning an NLK2pi cha
501 weight = bias_apply;
502 isElastic = false;
503 isStrangeProduction = true;
504 return new NNToNLK2piChannel(parti
505 } else if(NLKpiProductionCX>0.) {
506 INCL_WARN("Returning an NLKpi chan
507 weight = bias_apply;
508 isElastic = false;
509 isStrangeProduction = true;
510 return new NNToNLKpiChannel(partic
511 } else if(NLKProductionCX>0.) {
512 INCL_WARN("Returning an NLK channe
513 weight = bias_apply;
514 isElastic = false;
515 isStrangeProduction = true;
516 return new NNToNLKChannel(particle
517 } else if(omegafourPiProductionCX>0.)
518 INCL_WARN("Returning an Omega + fo
519 weight = counterweight;
520 isElastic = false;
521 return new NNOmegaToMultiPionsChan
522 } else if(omegathreePiProductionCX>0.)
523 INCL_WARN("Returning an Omega + th
524 weight = counterweight;
525 isElastic = false;
526 return new NNOmegaToMultiPionsChan
527 } else if(omegatwoPiProductionCX>0.) {
528 INCL_WARN("Returning an Omega + tw
529 weight = counterweight;
530 isElastic = false;
531 return new NNOmegaToMultiPionsChan
532 } else if(omegaonePiProductionCX>0.) {
533 INCL_WARN("Returning an Omega + on
534 weight = counterweight;
535 isElastic = false;
536 return new NNOmegaToMultiPionsChannel
537 } else if(omegadeltaProductionCX>0.) {
538 INCL_WARN("Returning an Omega + De
539 weight = counterweight;
540 isElastic = false;
541 return new NDeltaOmegaProductionCh
542 } else if(omegaProductionCX>0.) {
543 INCL_WARN("Returning an Omega chan
544 weight = counterweight;
545 isElastic = false;
546 return new NNToNNOmegaChannel(part
547 } else if(etafourPiProductionCX>0.) {
548 INCL_WARN("Returning an Eta + four
549 weight = counterweight;
550 isElastic = false;
551 return new NNEtaToMultiPionsChanne
552 } else if(etathreePiProductionCX>0.) {
553 INCL_WARN("Returning an Eta + thre
554 weight = counterweight;
555 isElastic = false;
556 return new NNEtaToMultiPionsChanne
557 } else if(etatwoPiProductionCX>0.) {
558 INCL_WARN("Returning an Eta + two
559 weight = counterweight;
560 isElastic = false;
561 return new NNEtaToMultiPionsChanne
562 } else if(etaonePiProductionCX>0.) {
563 INCL_WARN("Returning an Eta + one
564 weight = counterweight;
565 isElastic = false;
566 return new NNEtaToMultiPionsChanne
567 } else if(etadeltaProductionCX>0.) {
568 INCL_WARN("Returning an Eta + Delt
569 weight = counterweight;
570 isElastic = false;
571 return new NDeltaEtaProductionChan
572 } else if(etaProductionCX>0.) {
573 INCL_WARN("Returning an Eta channe
574 weight = counterweight;
575 isElastic = false;
576 return new NNToNNEtaChannel(partic
577 } else if(fourPiProductionCX>0.) {
578 INCL_WARN("Returning a 4pi channel
579 weight = counterweight;
580 isElastic = false;
581 return new NNToMultiPionsChannel(4
582 } else if(threePiProductionCX>0.) {
583 INCL_WARN("Returning a 3pi channel
584 weight = counterweight;
585 isElastic = false;
586 return new NNToMultiPionsChannel(3
587 } else if(twoPiProductionCX>0.) {
588 INCL_WARN("Returning a 2pi channel
589 weight = counterweight;
590 isElastic = false;
591 return new NNToMultiPionsChannel(2
592 } else if(onePiProductionCX>0.) {
593 INCL_WARN("Returning a 1pi channel
594 weight = counterweight;
595 isElastic = false;
596 return new NNToMultiPionsChannel(1
597 } else if(deltaProductionCX>0.) {
598 INCL_WARN("Returning a delta-produ
599 weight = counterweight;
600 isElastic = false;
601 return new DeltaProductionChannel(
602 } else {
603 INCL_WARN("Returning an elastic ch
604 weight = counterweight;
605 isElastic = true;
606 return new ElasticChannel(particle
607 }
608 }
609
610 //// NDelta
611 }
612 else if((particle1->isNucleon() && particl
613 (particle1->isDelta() && part
614
615 G4double NLKProductionCX = CrossSect
616 G4double NSKProductionCX = CrossSect
617 G4double DeltaLKProductionCX = Cross
618 G4double DeltaSKProductionCX = Cross
619 G4double NNKKbProductionCX = CrossSe
620
621 const G4double UnStrangeProdCX = Cro
622 const G4double StrangenessProdCX = (
623
624 G4double counterweight = (1. - bias_
625
626 if(counterweight < 0.5){
627 counterweight = 0.5;
628 bias_apply = 0.5*UnStrangeProdCX/
629
630 NLKProductionCX = CrossSections::
631 NSKProductionCX = CrossSections::
632 DeltaLKProductionCX = CrossSectio
633 DeltaSKProductionCX = CrossSectio
634 NNKKbProductionCX = CrossSections
635 }
636
637 G4double elasticCX = CrossSections::
638 G4double recombinationCX = CrossSect
639
640 const G4double rChannel=Random::shoo
641
642 if(elasticCX > rChannel) {
643 isElastic = true;
644 // Elastic N Delta channel
645 INCL_DEBUG("NDelta interaction: e
646 weight = counterweight;
647 return new ElasticChannel(particl
648 } else if (elasticCX + recombination
649 isElastic = false;
650 // Recombination
651 // NDelta -> NN channel is chosen
652 INCL_DEBUG("NDelta interaction: r
653 weight = counterweight;
654 return new RecombinationChannel(p
655 } else if (elasticCX + recombination
656 isElastic = false;
657 isStrangeProduction = true;
658 // NDelta -> NLK channel is chosen
659 INCL_DEBUG("NDelta interaction: N
660 weight = bias_apply;
661 return new NDeltaToNLKChannel(par
662 } else if (elasticCX + recombination
663 isElastic = false;
664 isStrangeProduction = true;
665 // NDelta -> NSK channel is chosen
666 INCL_DEBUG("NDelta interaction: N
667 weight = bias_apply;
668 return new NDeltaToNSKChannel(par
669 } else if (elasticCX + recombination
670 isElastic = false;
671 isStrangeProduction = true;
672 // NDelta -> DeltaLK channel is chosen
673 INCL_DEBUG("NDelta interaction: D
674 weight = bias_apply;
675 return new NDeltaToDeltaLKChannel
676 } else if (elasticCX + recombination
677 isElastic = false;
678 isStrangeProduction = true;
679 // NDelta -> DeltaSK channel is chosen
680 INCL_DEBUG("NDelta interaction: D
681 weight = bias_apply;
682 return new NDeltaToDeltaSKChannel
683 } else if (elasticCX + recombination
684 isElastic = false;
685 isStrangeProduction = true;
686 // NDelta -> NNKKb channel is chosen
687 INCL_DEBUG("NDelta interaction: N
688 weight = bias_apply;
689 return new NDeltaToNNKKbChannel(p
690 }
691 else{
692 INCL_ERROR("rChannel > (Strangene
693 weight = counterweight;
694 isElastic = true;
695 return new ElasticChannel(particl
696 }
697
698 //// DeltaDelta
699 } else if(particle1->isDelta() && particle
700 isElastic = true;
701 INCL_DEBUG("DeltaDelta interaction: el
702 return new ElasticChannel(particle1, p
703
704 //// PiN
705 } else if(isPiN) {
706
707 G4double LKProdCX = CrossSections::NpiTo
708 G4double SKProdCX = CrossSections::NpiTo
709 G4double LKpiProdCX = CrossSections::Npi
710 G4double SKpiProdCX = CrossSections::Npi
711 G4double LK2piProdCX = CrossSections::Np
712 G4double SK2piProdCX = CrossSections::Np
713 G4double NKKbProdCX = CrossSections::Npi
714 G4double MissingCX = CrossSections::NpiT
715
716 const G4double UnStrangeProdCX = CrossSe
717 + CrossSect
718 + CrossSect
719 const G4double StrangenessProdCX = (LKPr
720
721 G4double counterweight = (1. - bias_appl
722
723 if(counterweight < 0.5) {
724 counterweight = 0.5;
725 bias_apply = 0.5*UnStrangeProdCX/Stra
726 LKProdCX = CrossSections::NpiToLK(par
727 SKProdCX = CrossSections::NpiToSK(par
728 LKpiProdCX = CrossSections::NpiToLKpi
729 SKpiProdCX = CrossSections::NpiToSKpi
730 LK2piProdCX = CrossSections::NpiToLK2
731 SK2piProdCX = CrossSections::NpiToSK2
732 NKKbProdCX = CrossSections::NpiToNKKb
733 MissingCX = CrossSections::NpiToMissi
734 }
735
736
737 const G4double elasticCX = CrossSections
738 const G4double deltaProductionCX = Cross
739 const G4double onePiProductionCX = Cross
740 const G4double twoPiProductionCX = Cross
741 const G4double threePiProductionCX =
742 const G4double etaProductionCX = Cros
743 const G4double omegaProductionCX = Cr
744
745 const G4double totCX=CrossSections::tota
746
747 // assert(std::fabs(totCX-elasticCX-deltaProdu
748
749 const G4double rChannel=Random::shoot()
750
751 if(elasticCX > rChannel) {
752 isElastic = true;
753 // Elastic PiN channel
754 INCL_DEBUG("PiN interaction: elastic c
755 weight = counterweight;
756 return new PiNElasticChannel(particle1
757 } else if(elasticCX + deltaProductionCX
758 isElastic = false;
759 // PiN -> Delta channel is chosen
760 INCL_DEBUG("PiN interaction: Delta cha
761 weight = counterweight;
762 return new PiNToDeltaChannel(particle1
763 } else if(elasticCX + deltaProductionCX
764 isElastic = false;
765 // PiN -> PiNPi channel is chosen
766 INCL_DEBUG("PiN interaction: one Pion
767 weight = counterweight;
768 return new PiNToMultiPionsChannel(2,pa
769 } else if(elasticCX + deltaProductionCX
770 isElastic = false;
771 // PiN -> PiN2Pi channel is chosen
772 INCL_DEBUG("PiN interaction: two Pions
773 weight = counterweight;
774 return new PiNToMultiPionsChannel(3,pa
775 } else if(elasticCX + deltaProductionCX
776 isElastic = false;
777 // PiN -> PiN3Pi channel is chosen
778 INCL_DEBUG("PiN interaction: three Pio
779 weight = counterweight;
780 return new PiNToMultiPionsChannel(4,pa
781 } else if(elasticCX + deltaProductionCX
782 isElastic = false;
783 // PiN -> EtaN channel is chosen
784 INCL_DEBUG("PiN interaction: Eta chann
785 weight = counterweight;
786 return new PiNToEtaChannel(particle1,
787 } else if(elasticCX + deltaProductionCX
788 isElastic = false;
789 // PiN -> OmegaN channel is chosen
790 INCL_DEBUG("PiN interaction: Omega cha
791 weight = counterweight;
792 return new PiNToOmegaChannel(particle1
793 } else if(elasticCX + deltaProductionCX
794 + LKProdCX > rChanne
795 isElastic = false;
796 isStrangeProduction = true;
797 // PiN -> LK channel is chosen
798 INCL_DEBUG("PiN interaction: LK channe
799 weight = bias_apply;
800 return new NpiToLKChannel(particle1, p
801 } else if(elasticCX + deltaProductionCX
802 + LKProdCX + SKProdC
803 isElastic = false;
804 isStrangeProduction = true;
805 // PiN -> SK channel is chosen
806 INCL_DEBUG("PiN interaction: SK channe
807 weight = bias_apply;
808 return new NpiToSKChannel(particle1, p
809 } else if(elasticCX + deltaProductionCX
810 + LKProdCX + SKProdC
811 isElastic = false;
812 isStrangeProduction = true;
813 // PiN -> LKpi channel is chosen
814 INCL_DEBUG("PiN interaction: LKpi chan
815 weight = bias_apply;
816 return new NpiToLKpiChannel(particle1,
817 } else if(elasticCX + deltaProductionCX
818 + LKProdCX + SKProdC
819 isElastic = false;
820 isStrangeProduction = true;
821 // PiN -> SKpi channel is chosen
822 INCL_DEBUG("PiN interaction: SKpi chan
823 weight = bias_apply;
824 return new NpiToSKpiChannel(particle1,
825 } else if(elasticCX + deltaProductionCX
826 + LKProdCX + SKProdC
827 isElastic = false;
828 isStrangeProduction = true;
829 // PiN -> LK2pi channel is chosen
830 INCL_DEBUG("PiN interaction: LK2pi cha
831 weight = bias_apply;
832 return new NpiToLK2piChannel(particle1
833 } else if(elasticCX + deltaProductionCX
834 + LKProdCX + SKProdC
835 isElastic = false;
836 isStrangeProduction = true;
837 // PiN -> SK2pi channel is chosen
838 INCL_DEBUG("PiN interaction: SK2pi cha
839 weight = bias_apply;
840 return new NpiToSK2piChannel(particle1
841 } else if(elasticCX + deltaProductionCX
842 + LKProdCX + SKProdC
843 isElastic = false;
844 isStrangeProduction = true;
845 // PiN -> NKKb channel is chosen
846 INCL_DEBUG("PiN interaction: NKKb chan
847 weight = bias_apply;
848 return new NpiToNKKbChannel(particle1,
849 } else if(elasticCX + deltaProductionCX
850 + LKProdCX + SKProdC
851 isElastic = false;
852 isStrangeProduction = true;
853 // PiN -> Missinge Strangeness channel is chos
854 INCL_DEBUG("PiN interaction: Missinge
855 weight = bias_apply;
856 return new NpiToMissingStrangenessChan
857 }
858 else {
859 INCL_WARN("inconsistency within the P
860 if(MissingCX>0.) {
861 INCL_WARN("Returning a Missinge St
862 weight = bias_apply;
863 isElastic = false;
864 isStrangeProduction = true;
865 return new NpiToMissingStrangeness
866 } else if(NKKbProdCX>0.) {
867 INCL_WARN("Returning a NKKb channe
868 weight = bias_apply;
869 isElastic = false;
870 isStrangeProduction = true;
871 return new NpiToNKKbChannel(partic
872 } else if(SK2piProdCX>0.) {
873 INCL_WARN("Returning a SK2pi chann
874 weight = bias_apply;
875 isElastic = false;
876 isStrangeProduction = true;
877 return new NpiToSK2piChannel(parti
878 } else if(LK2piProdCX>0.) {
879 INCL_WARN("Returning a LK2pi chann
880 weight = bias_apply;
881 isElastic = false;
882 isStrangeProduction = true;
883 return new NpiToLK2piChannel(parti
884 } else if(SKpiProdCX>0.) {
885 INCL_WARN("Returning a SKpi channe
886 weight = bias_apply;
887 isElastic = false;
888 isStrangeProduction = true;
889 return new NpiToSKpiChannel(partic
890 } else if(LKpiProdCX>0.) {
891 INCL_WARN("Returning a LKpi channe
892 weight = bias_apply;
893 isElastic = false;
894 isStrangeProduction = true;
895 return new NpiToLKpiChannel(partic
896 } else if(SKProdCX>0.) {
897 INCL_WARN("Returning a SK channel"
898 weight = bias_apply;
899 isElastic = false;
900 isStrangeProduction = true;
901 return new NpiToSKChannel(particle
902 } else if(LKProdCX>0.) {
903 INCL_WARN("Returning a LK channel"
904 weight = bias_apply;
905 isElastic = false;
906 isStrangeProduction = true;
907 return new NpiToLKChannel(particle
908 } else if(omegaProductionCX>0.) {
909 INCL_WARN("Returning a Omega chann
910 weight = counterweight;
911 isElastic = false;
912 return new PiNToOmegaChannel(parti
913 } else if(etaProductionCX>0.) {
914 INCL_WARN("Returning a Eta channel
915 weight = counterweight;
916 isElastic = false;
917 return new PiNToEtaChannel(particl
918 } else if(threePiProductionCX>0.) {
919 INCL_WARN("Returning a 3pi channel
920 weight = counterweight;
921 isElastic = false;
922 return new PiNToMultiPionsChannel(
923 } else if(twoPiProductionCX>0.) {
924 INCL_WARN("Returning a 2pi channel
925 weight = counterweight;
926 isElastic = false;
927 return new PiNToMultiPionsChannel(
928 } else if(onePiProductionCX>0.) {
929 INCL_WARN("Returning a 1pi channel
930 weight = counterweight;
931 isElastic = false;
932 return new PiNToMultiPionsChannel(
933 } else if(deltaProductionCX>0.) {
934 INCL_WARN("Returning a delta-produ
935 weight = counterweight;
936 isElastic = false;
937 return new PiNToDeltaChannel(parti
938 } else {
939 INCL_WARN("Returning an elastic ch
940 weight = counterweight;
941 isElastic = true;
942 return new PiNElasticChannel(parti
943 }
944 }
945 } else if ((particle1->isNucleon() && part
946 //// EtaN
947
948 const G4double elasticCX =
949 const G4double onePiProduc
950 const G4double twoPiProduc
951 const G4double totCX=Cross
952 // assert(std::fabs(totCX-elasticCX-onePiProdu
953
954 const G4double rChannel=Ra
955
956 if(elasticCX > rChannel) {
957 // Elastic EtaN channel
958 isElastic = true;
959 INCL_DEBUG("EtaN inter
960 return new EtaNElastic
961 } else if(elasticCX + oneP
962 isElastic = false;
963 // EtaN -> EtaPiN channel is chosen
964 INCL_DEBUG("EtaN inter
965 return new EtaNToPiNCh
966 } else if(elasticCX + oneP
967 isElastic = false;
968 // EtaN -> EtaPiPiN channel is chosen
969 INCL_DEBUG("EtaN inter
970 return new EtaNToPiPiN
971 }
972
973 else {
974 INCL_WARN("inconsisten
975 if(twoPiProductionCX>0
976 INCL_WARN("Returni
977 isElastic = false;
978 return new EtaNToP
979 } else if(onePiProduct
980 INCL_WARN("Returni
981 isElastic = false;
982 return new EtaNToP
983 } else {
984 INCL_WARN("Returni
985 isElastic = true;
986 return new EtaNEla
987 }
988 }
989
990 } else if ((particle1->isNucleon() && part
991 //// OmegaN
992
993 const G4double elasticCX =
994 const G4double onePiProduc
995 const G4double twoPiProduc
996 const G4double totCX=Cross
997 // assert(std::fabs(totCX-elasticCX-onePiProdu
998
999 const G4double rChannel=Ra
1000
1001 if(elasticCX > rChannel)
1002 // Elastic OmegaN channel
1003 isElastic = true;
1004 INCL_DEBUG("OmegaN in
1005 return new OmegaNElas
1006 } else if(elasticCX + one
1007 isElastic = false;
1008 // OmegaN -> PiN channel is chosen
1009 INCL_DEBUG("OmegaN interaction: P
1010 return new OmegaNToPiNChannel(par
1011 } else if(elasticCX + one
1012 isElastic = false;
1013 // OmegaN -> PiPiN channel is chosen
1014 INCL_DEBUG("OmegaN in
1015 return new OmegaNToPi
1016 }
1017 else {
1018 INCL_WARN("inconsiste
1019 if(twoPiProductionCX>
1020 INCL_WARN("Return
1021 isElastic = false
1022 return new OmegaN
1023 } else if(onePiProduc
1024 INCL_WARN("Return
1025 isElastic = false
1026 return new OmegaN
1027 } else {
1028 INCL_WARN("Return
1029 isElastic = true;
1030 return new OmegaN
1031 }
1032 }
1033 } else if ((particle1->isNucleon() && par
1034 //// KN
1035 const G4double elasticCX = CrossSecti
1036 const G4double quasielasticCX = Cross
1037 const G4double NKToNKpiCX = CrossSect
1038 const G4double NKToNK2piCX = CrossSec
1039 const G4double totCX=CrossSections::t
1040 // assert(std::fabs(totCX-elasticCX-quasielas
1041
1042 const G4double rChannel=Random::shoot
1043 if(elasticCX > rChannel){
1044 // Elastic KN channel is chosen
1045 isElastic = true;
1046 INCL_DEBUG("KN interaction: elast
1047 return new NKElasticChannel(parti
1048 } else if(elasticCX + quasielasticCX
1049 // Quasi-elastic KN channel is chosen
1050 isElastic = false; // true ??
1051 INCL_DEBUG("KN interaction: quasi
1052 return new NKToNKChannel(particle
1053 } else if(elasticCX + quasielasticCX
1054 // KN -> NKpi channel is chosen
1055 isElastic = false;
1056 INCL_DEBUG("KN interaction: NKpi
1057 return new NKToNKpiChannel(partic
1058 } else if(elasticCX + quasielasticCX
1059 // KN -> NK2pi channel is chosen
1060 isElastic = false;
1061 INCL_DEBUG("KN interaction: NK2pi
1062 return new NKToNK2piChannel(parti
1063 } else {
1064 INCL_WARN("inconsistency within t
1065 if(NKToNK2piCX>0.) {
1066 INCL_WARN("Returning a NKToNK
1067 isElastic = false;
1068 return new NKToNK2piChannel(p
1069 } else if(NKToNKpiCX>0.) {
1070 INCL_WARN("Returning a NKToNK
1071 isElastic = false;
1072 return new NKToNKpiChannel(pa
1073 } else if(quasielasticCX>0.) {
1074 INCL_WARN("Returning a quasi-
1075 isElastic = false; // true ??
1076 return new NKToNKChannel(part
1077 } else {
1078 INCL_WARN("Returning an elast
1079 isElastic = true;
1080 return new NKElasticChannel(p
1081 }
1082 }
1083 } else if ((particle1->isNucleon() && par
1084 //// KbN
1085 const G4double elasticCX = CrossSecti
1086 const G4double quasielasticCX = Cross
1087 const G4double NKbToNKbpiCX = CrossSe
1088 const G4double NKbToNKb2piCX = CrossS
1089 const G4double NKbToLpiCX = CrossSect
1090 const G4double NKbToL2piCX = CrossSec
1091 const G4double NKbToSpiCX = CrossSect
1092 const G4double NKbToS2piCX = CrossSec
1093 const G4double totCX=CrossSections::t
1094 // assert(std::fabs(totCX-elasticCX-quasielas
1095
1096 const G4double rChannel=Random::shoot
1097 if(elasticCX > rChannel){
1098 // Elastic KbN channel is chosen
1099 isElastic = true;
1100 INCL_DEBUG("KbN interaction: elas
1101 return new NKbElasticChannel(part
1102 } else if(elasticCX + quasielasticCX
1103 // Quasi-elastic KbN channel is chosen
1104 isElastic = false; // true ??
1105 INCL_DEBUG("KbN interaction: quas
1106 return new NKbToNKbChannel(partic
1107 } else if(elasticCX + quasielasticCX
1108 // KbN -> NKbpi channel is chosen
1109 isElastic = false;
1110 INCL_DEBUG("KbN interaction: NKbp
1111 return new NKbToNKbpiChannel(part
1112 } else if(elasticCX + quasielasticCX
1113 // KbN -> NKb2pi channel is chosen
1114 isElastic = false;
1115 INCL_DEBUG("KbN interaction: NKb2
1116 return new NKbToNKb2piChannel(par
1117 } else if(elasticCX + quasielasticCX
1118 // KbN -> Lpi channel is chosen
1119 isElastic = false;
1120 INCL_DEBUG("KbN interaction: Lpi
1121 return new NKbToLpiChannel(partic
1122 } else if(elasticCX + quasielasticCX
1123 // KbN -> L2pi channel is chosen
1124 isElastic = false;
1125 INCL_DEBUG("KbN interaction: L2pi
1126 return new NKbToL2piChannel(parti
1127 } else if(elasticCX + quasielasticCX
1128 // KbN -> Spi channel is chosen
1129 isElastic = false;
1130 INCL_DEBUG("KbN interaction: Spi
1131 return new NKbToSpiChannel(partic
1132 } else if(elasticCX + quasielasticCX
1133 // KbN -> S2pi channel is chosen
1134 isElastic = false;
1135 INCL_DEBUG("KbN interaction: S2pi
1136 return new NKbToS2piChannel(parti
1137 } else {
1138 INCL_WARN("inconsistency within t
1139 if(NKbToS2piCX>0.) {
1140 INCL_WARN("Returning a NKbToS
1141 isElastic = false;
1142 return new NKbToS2piChannel(p
1143 } else if(NKbToSpiCX>0.) {
1144 INCL_WARN("Returning a NKbToS
1145 isElastic = false;
1146 return new NKbToSpiChannel(pa
1147 } else if(NKbToL2piCX>0.) {
1148 INCL_WARN("Returning a NKbToL
1149 isElastic = false;
1150 return new NKbToL2piChannel(p
1151 } else if(NKbToLpiCX>0.) {
1152 INCL_WARN("Returning a NKbToL
1153 isElastic = false;
1154 return new NKbToLpiChannel(pa
1155 } else if(NKbToNKb2piCX>0.) {
1156 INCL_WARN("Returning a NKbToN
1157 isElastic = false;
1158 return new NKbToNKb2piChannel
1159 } else if(NKbToNKbpiCX>0.) {
1160 INCL_WARN("Returning a NKbToN
1161 isElastic = false;
1162 return new NKbToNKbpiChannel(
1163 } else if(quasielasticCX>0.) {
1164 INCL_WARN("Returning a quasi-
1165 isElastic = false; // true ??
1166 return new NKbToNKbChannel(pa
1167 } else {
1168 INCL_WARN("Returning an elast
1169 isElastic = true;
1170 return new NKbElasticChannel(
1171 }
1172 }
1173 } else if ((particle1->isNucleon() && par
1174 //// NLambda
1175 const G4double elasticCX = CrossSecti
1176 const G4double NLToNSCX = CrossSectio
1177 const G4double totCX=CrossSections::t
1178 // assert(std::fabs(totCX-elasticCX-NLToNSCX)
1179
1180 const G4double rChannel=Random::shoot
1181 if(elasticCX > rChannel){
1182 // Elastic NLambda channel is chosen
1183 isElastic = true;
1184 INCL_DEBUG("NLambda interaction:
1185 return new NYElasticChannel(parti
1186 } else if(elasticCX + NLToNSCX > rCha
1187 // Quasi-elastic NLambda channel is chosen
1188 isElastic = false; // true ??
1189 INCL_DEBUG("NLambda interaction:
1190 return new NLToNSChannel(particle
1191 } else {
1192 INCL_WARN("inconsistency within t
1193 if(NLToNSCX>0.) {
1194 INCL_WARN("Returning a quasi-
1195 isElastic = false; // true ??
1196 return new NLToNSChannel(part
1197 } else {
1198 INCL_WARN("Returning an elast
1199 isElastic = true;
1200 return new NYElasticChannel(p
1201 }
1202 }
1203 } else if ((particle1->isNucleon() && par
1204 //// NSigma
1205 const G4double elasticCX = CrossSecti
1206 const G4double NSToNLCX = CrossSectio
1207 const G4double NSToNSCX = CrossSectio
1208 const G4double totCX=CrossSections::t
1209 // assert(std::fabs(totCX-elasticCX-NSToNLCX-
1210
1211 const G4double rChannel=Random::shoot
1212 if(elasticCX > rChannel){
1213 // Elastic NSigma channel is chosen
1214 isElastic = true;
1215 INCL_DEBUG("NSigma interaction: e
1216 return new NYElasticChannel(parti
1217 } else if(elasticCX + NSToNLCX > rCha
1218 // NSigma -> NLambda channel is chosen
1219 isElastic = false; // true ??
1220 INCL_DEBUG("NSigma interaction: N
1221 return new NSToNLChannel(particle
1222 } else if(elasticCX + NSToNLCX + NSTo
1223 // NSigma -> NSigma quasi-elastic channel is
1224 isElastic = false; // true ??
1225 INCL_DEBUG("NSigma interaction: N
1226 return new NSToNSChannel(particle
1227 } else {
1228 INCL_WARN("inconsistency within t
1229 if(NSToNSCX>0.) {
1230 INCL_WARN("Returning a quasi-
1231 isElastic = false; // true ??
1232 return new NSToNSChannel(part
1233 } else if(NSToNLCX>0.) {
1234 INCL_WARN("Returning a NLambd
1235 isElastic = false; // true ??
1236 return new NSToNLChannel(part
1237 } else {
1238 INCL_WARN("Returning an elast
1239 isElastic = true;
1240 return new NYElasticChannel(p
1241 }
1242 }
1243 } else if ((particle1->isNucleon() && par
1244 //// NNbar
1245 const G4double totCX = CrossSections:
1246 const G4double NNbElasticCX = CrossSe
1247 const G4double NNbCEXCX = CrossSectio
1248 const G4double NNbToLLbCX = CrossSect
1249 const G4double NNbToNNbpiCX = CrossSe
1250 const G4double NNbToNNb2piCX = CrossS
1251 const G4double NNbToNNb3piCX = CrossS
1252 const G4double AnnihilationCX = Cross
1253
1254 // assert(std::fabs(totCX-NNbElasticCX-NNbCEX
1255
1256 const G4double rChannel=Random::shoot
1257 if (NNbElasticCX > rChannel) {
1258 // NNbar (elastic) channel is cho
1259 isElastic = true;
1260 //INCL_WARN("NNbar interaction: N
1261 return new NNbarElasticChannel(pa
1262 } else if (NNbElasticCX + NNbCEXCX >
1263 // NNbar (CEX) channel is chosen
1264 isElastic = false; // may be char
1265 //INCL_WARN("NNbar interaction: N
1266 return new NNbarCEXChannel(partic
1267 } else if (NNbElasticCX + NNbCEXCX +
1268 // NNbarToLLbar channel is chosen
1269 isElastic = false; // may be char
1270 //INCL_WARN("NNbar interaction: N
1271 return new NNbarToLLbarChannel(pa
1272 } else if (NNbElasticCX + NNbCEXCX +
1273 // NNbar to NNbar pi channel is c
1274 isElastic = false;
1275 //INCL_WARN("NNbar interaction: N
1276 return new NNbarToNNbarpiChannel(
1277 } else if (NNbElasticCX + NNbCEXCX +
1278 // NNbar to NNbar 2pi channel is
1279 isElastic = false;
1280 //INCL_WARN("NNbar interaction: N
1281 return new NNbarToNNbar2piChannel
1282 } else if (NNbElasticCX + NNbCEXCX +
1283 // NNbar to NNbar 3pi channel is
1284 isElastic = false;
1285 //INCL_WARN("NNbar interaction: N
1286 return new NNbarToNNbar3piChannel
1287 } else if (NNbElasticCX + NNbCEXCX +
1288 // NNbar annihilation channel is
1289 isElastic = false;
1290 AnnihilationType atype;
1291 if((particle1->getType()==antiPro
1292 atype = PTypeInFlight;
1293 }
1294 else if((particle1->getType()==an
1295 atype = NTypeInFlight;
1296 }
1297 else if((particle1->getType()==an
1298 atype = NbarPTypeInFlight;
1299 }
1300 else if((particle1->getType()==an
1301 atype = NbarNTypeInFlight;
1302 }
1303 else{
1304 atype = Def;
1305 INCL_ERROR("Annihilation type p
1306 }
1307 theNucleus->setAType(atype);
1308 return new NNbarToAnnihilationCha
1309 } else {
1310 INCL_WARN("Inconsistency within
1311 if (NNbToNNb3piCX > 0.0) {
1312 INCL_WARN("Returning an NNb
1313 isElastic = false;
1314 return new NNbarToNNbar3piC
1315 } else if (NNbToNNb2piCX > 0.0)
1316 INCL_WARN("Returning an NNb
1317 isElastic = false;
1318 return new NNbarToNNbar2piC
1319 } else if (NNbToNNbpiCX > 0.0)
1320 INCL_WARN("Returning an NNb
1321 isElastic = false;
1322 return new NNbarToNNbarpiCh
1323 } else if (AnnihilationCX > 0.0
1324 INCL_WARN("Returning an NNb
1325 isElastic = false;
1326 AnnihilationType atype;
1327 if((particle1->getType()==a
1328 atype = PTypeInFlight;
1329 }
1330 else if((particle1->getType
1331 atype = NTypeInFlight;
1332 }
1333 else if((particle1->getType
1334 atype = NbarPTypeInFlight
1335 }
1336 else if((particle1->getType
1337 atype = NbarNTypeInFlight
1338 }
1339 else{
1340 atype = Def;
1341 INCL_ERROR("Annihilation
1342 }
1343 theNucleus->setAType(atype)
1344 return new NNbarToAnnihilat
1345 } else if (NNbCEXCX > 0.0) {
1346 INCL_WARN("Returning an NNb
1347 isElastic = false;
1348 return new NNbarCEXChannel(
1349 } else if (NNbToLLbCX > 0.0) {
1350 INCL_WARN("Returning an NNb
1351 isElastic = false;
1352 return new NNbarToLLbarChan
1353 } else {
1354 INCL_WARN("Elastic NNbar ch
1355 isElastic = true;
1356 return new NNbarElasticChan
1357 }
1358 }
1359 }
1360
1361 else {
1362 INCL_DEBUG("BinaryCollisionAvatar can o
1363 << '\n'
1364 << particle1->print()
1365 << '\n'
1366 << particle2->print()
1367 << '\n');
1368 InteractionAvatar::restoreParticles();
1369 return NULL;
1370 }
1371 }
1372
1373 void BinaryCollisionAvatar::preInteraction(
1374 isParticle1Spectator = particle1->isTarge
1375 isParticle2Spectator = particle2->isTarge
1376 InteractionAvatar::preInteraction();
1377 }
1378
1379 void BinaryCollisionAvatar::postInteraction
1380 // Call the postInteraction method of the
1381 // (provides Pauli blocking and enforces
1382 InteractionAvatar::postInteraction(fs);
1383
1384 switch(fs->getValidity()) {
1385 case PauliBlockedFS:
1386 theNucleus->getStore()->getBook().inc
1387 break;
1388 case NoEnergyConservationFS:
1389 case ParticleBelowFermiFS:
1390 case ParticleBelowZeroFS:
1391 break;
1392 case ValidFS:
1393 Book &theBook = theNucleus->getStore(
1394 theBook.incrementAcceptedCollisions()
1395
1396 if(theBook.getAcceptedCollisions() ==
1397 // Store time and cross section of
1398 G4double t = theBook.getCurrentTime
1399 theBook.setFirstCollisionTime(t);
1400 theBook.setFirstCollisionXSec(oldXS
1401 // Increase the number of Kaon by 1
1402 if(isStrangeProduction) theNucleus-
1403 // Store position and momentum of t
1404 // collision
1405 if((isParticle1Spectator && isParti
1406 INCL_ERROR("First collision must
1407 }
1408 if(isParticle1Spectator) {
1409 theBook.setFirstCollisionSpectato
1410 theBook.setFirstCollisionSpectato
1411 } else {
1412 theBook.setFirstCollisionSpectato
1413 theBook.setFirstCollisionSpectato
1414 }
1415
1416 // Store the elasticity of the firs
1417 theBook.setFirstCollisionIsElastic(
1418 }
1419 }
1420 return;
1421 }
1422
1423 std::string BinaryCollisionAvatar::dump() c
1424 std::stringstream ss;
1425 ss << "(avatar " << theTime <<" 'nn-colli
1426 << "(list " << '\n'
1427 << particle1->dump()
1428 << particle2->dump()
1429 << "))" << '\n';
1430 return ss.str();
1431 }
1432
1433 }
1434