Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCL::BinaryCollisionAvatar Class Reference

#include <G4INCLBinaryCollisionAvatar.hh>

Inheritance diagram for G4INCL::BinaryCollisionAvatar:

Public Member Functions

 BinaryCollisionAvatar (G4double, G4double, G4INCL::Nucleus *, G4INCL::Particle *, G4INCL::Particle *)
virtual ~BinaryCollisionAvatar ()
G4INCL::IChannel * getChannel ()
ParticleList getParticles () const
virtual void preInteraction ()
virtual void postInteraction (FinalState *)
std::string dump () const
Public Member Functions inherited from G4INCL::InteractionAvatar
 InteractionAvatar (G4double, G4INCL::Nucleus *, G4INCL::Particle *)
 InteractionAvatar (G4double, G4INCL::Nucleus *, G4INCL::Particle *, G4INCL::Particle *)
virtual ~InteractionAvatar ()
void setSrcPartner (Particle *p)
void preInteractionLocalEnergy (Particle *const p)
 Apply local-energy transformation, if appropriate.
ThreeVector getboostVector ()
void setboostVector (ThreeVector &v)
Public Member Functions inherited from G4INCL::IAvatar
 IAvatar ()
 IAvatar (G4double time)
virtual ~IAvatar ()
FinalState * getFinalState ()
void fillFinalState (FinalState *fs)
G4double getTime () const
G4double setTime (G4double t) const
AvatarType getType () const
G4bool isACollision () const
G4bool isADecay () const
void setType (AvatarType t)
long getID () const
std::string toString ()

Static Public Member Functions

static void setCutNN (const G4double c)
static G4double getCutNN ()
static G4double getCutNNSquared ()
static G4double getBias ()
 Get the global bias factor.
static void setBias (const G4double b)
 Set the global bias factor.
Static Public Member Functions inherited from G4INCL::InteractionAvatar
static void deleteBackupParticles ()
 Release the memory allocated for the backup particles.
static InteractionAvatar * Instance ()

Additional Inherited Members

Static Public Attributes inherited from G4INCL::InteractionAvatar
static const G4double locEAccuracy = 1.E-4
 Target accuracy in the determination of the local-energy Q-value.
static const G4int maxIterLocE = 50
 Max number of iterations for the determination of the local-energy Q-value.
Protected Member Functions inherited from G4INCL::InteractionAvatar
G4bool bringParticleInside (Particle *const p)
void preInteractionBlocking ()
 Store the state of the particles before the interaction.
void preInteraction ()
void postInteraction (FinalState *)
void restoreParticles () const
 Restore the state of both particles.
void restoreSrcPartner (FinalState *fs)
G4bool shouldUseLocalEnergy () const
 true if the given avatar should use local energy
G4bool enforceEnergyConservation (FinalState *const fs)
 Enforce energy conservation.
Protected Attributes inherited from G4INCL::InteractionAvatar
EventInfo theEventInfo
Nucleus * theNucleus
Particle * particle1
Particle * particle2
ThreeVector boostVector
G4double oldTotalEnergy
G4double oldXSec
G4bool isPiN
G4double weight
ParticleList modified
ParticleList created
ParticleList modifiedAndCreated
ParticleList Destroyed
ParticleList ModifiedAndDestroyed
Protected Attributes inherited from G4INCL::IAvatar
G4double theTime
Static Protected Attributes inherited from G4INCL::InteractionAvatar
static G4ThreadLocal Particle * backupParticle1 = NULL
static G4ThreadLocal Particle * backupParticle2 = NULL

Detailed Description

Definition at line 57 of file G4INCLBinaryCollisionAvatar.hh.

Constructor & Destructor Documentation

◆ BinaryCollisionAvatar()

G4INCL::BinaryCollisionAvatar::BinaryCollisionAvatar ( G4double time,
G4double crossSection,
G4INCL::Nucleus * n,
G4INCL::Particle * p1,
G4INCL::Particle * p2 )

Definition at line 133 of file G4INCLBinaryCollisionAvatar.cc.

135 : InteractionAvatar(time, n, p1, p2), theCrossSection(crossSection),
136 isParticle1Spectator(false),
137 isParticle2Spectator(false),
138 isElastic(false),
139 isStrangeProduction(false)
140 {
141 setType(CollisionAvatarType);
142 }
G4INCL::IAvatar::setType
void setType(AvatarType t)
Definition G4INCLIAvatar.hh:94
G4INCL::InteractionAvatar::InteractionAvatar
InteractionAvatar(G4double, G4INCL::Nucleus *, G4INCL::Particle *)
Definition G4INCLInteractionAvatar.cc:71
G4INCL::CollisionAvatarType
@ CollisionAvatarType
Definition G4INCLIAvatar.hh:67

◆ ~BinaryCollisionAvatar()

G4INCL::BinaryCollisionAvatar::~BinaryCollisionAvatar ( )
virtual

Definition at line 144 of file G4INCLBinaryCollisionAvatar.cc.

144 {
145 }

Member Function Documentation

◆ dump()

std::string G4INCL::BinaryCollisionAvatar::dump ( ) const
virtual

Implements G4INCL::IAvatar.

Definition at line 1630 of file G4INCLBinaryCollisionAvatar.cc.

1630 {
1631 std::stringstream ss;
1632 ss << "(avatar " << theTime <<" 'nn-collision" << '\n'
1633 << "(list " << '\n'
1634 << particle1->dump()
1635 << particle2->dump()
1636 << "))" << '\n';
1637 return ss.str();
1638 }

◆ getBias()

G4double G4INCL::BinaryCollisionAvatar::getBias ( )
inlinestatic

Get the global bias factor.

Definition at line 84 of file G4INCLBinaryCollisionAvatar.hh.

84{ return bias; }

◆ getChannel()

G4INCL::IChannel * G4INCL::BinaryCollisionAvatar::getChannel ( )
virtual

Check again the distance of approach. In order for the avatar to be realised, we have to perform a check in the CM system. We define a distance four-vector as

\‍[ (0, \Delta\vec{x}), \‍]

where $\Delta\vec{x}$ is the distance vector of the particles at their minimum distance of approach (i.e. at the avatar time). By boosting this four-vector to the CM frame of the two particles and we obtain a new four vector

\‍[ (\Delta t', \Delta\vec{x}'), \‍]

with a non-zero time component (the collision happens simultaneously for the two particles in the lab system, but not in the CM system). In order for the avatar to be realised, we require that

\‍[ |\Delta\vec{x}'| \leq \sqrt{\sigma/\pi}.\‍]

Note that $|\Delta\vec{x}'|\leq|\Delta\vec{x}|$; thus, the condition above is more restrictive than the check that we perform in G4INCL::Propagation::StandardPropagationModel::generateBinaryCollisionAvatar. In other words, the avatar generation cannot miss any physical collision avatars.

Bias apply for this reaction in order to get the same ParticleBias for all stange particles. Can be reduced after because of the safeguard.

Implements G4INCL::InteractionAvatar.

Definition at line 208 of file G4INCLBinaryCollisionAvatar.cc.

208 {
209 // We already check cutNN at avatar creation time, but we have to check it
210 // again here. For composite projectiles, we might have created independent
211 // avatars with no cutNN before any collision took place.
212 if(particle1->isNucleon()
213 && particle2->isNucleon()
214 && theNucleus->getStore()->getBook().getAcceptedCollisions()!=0) {
215 const G4double energyCM2 = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
216 // Below a certain cut value we don't do anything:
217 if(energyCM2 < cutNNSquared) {
218 INCL_DEBUG("CM energy = sqrt(" << energyCM2 << ") MeV < std::sqrt(" << cutNNSquared
219 << ") MeV = cutNN" << "; returning a NULL channel" << '\n');
220 InteractionAvatar::restoreParticles();
221 return NULL;
222 }
223 }
224
225 /** Check again the distance of approach. In order for the avatar to be
226 * realised, we have to perform a check in the CM system. We define a
227 * distance four-vector as
228 * \f[ (0, \Delta\vec{x}), \f]
229 * where \f$\Delta\vec{x}\f$ is the distance vector of the particles at
230 * their minimum distance of approach (i.e. at the avatar time). By
231 * boosting this four-vector to the CM frame of the two particles and we
232 * obtain a new four vector
233 * \f[ (\Delta t', \Delta\vec{x}'), \f]
234 * with a non-zero time component (the collision happens simultaneously for
235 * the two particles in the lab system, but not in the CM system). In order
236 * for the avatar to be realised, we require that
237 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/\pi}.\f]
238 * Note that \f$|\Delta\vec{x}'|\leq|\Delta\vec{x}|\f$; thus, the condition
239 * above is more restrictive than the check that we perform in
240 * G4INCL::Propagation::StandardPropagationModel::generateBinaryCollisionAvatar.
241 * In other words, the avatar generation cannot miss any physical collision
242 * avatars.
243 */
244 ThreeVector minimumDistance = particle1->getPosition();
245 minimumDistance -= particle2->getPosition();
246 const G4double betaDotX = boostVector.dot(minimumDistance);
247 const G4double minDist = Math::tenPi*(minimumDistance.mag2() + betaDotX*betaDotX / (1.-boostVector.mag2()));
248
249 Config const *theConfig=theNucleus->getStore()->getConfig();
250 if ((minDist > theCrossSection) &&
251 (!((particle1->getType()==antiProton && particle1->getKineticEnergy() <= theConfig->getAtrestThreshold()) ||
252 (particle2->getType()==antiProton && particle2->getKineticEnergy() <= theConfig->getAtrestThreshold()) ||
253 (particle1->getType()==antiNeutron && particle1->getKineticEnergy() <= particle1->getPotentialEnergy()) ||
254 (particle2->getType()==antiNeutron && particle2->getKineticEnergy() <= particle2->getPotentialEnergy())))) {
255 if(!((particle1->isAntiNucleon() && particle1->getEnergy() <= particle1->getINCLMass()) ||
256 (particle2->isAntiNucleon() && particle2->getEnergy() <= particle2->getINCLMass()))){
257 INCL_DEBUG("CM distance of approach is too small: " << minDist << ">" <<
258 theCrossSection <<"; returning a NULL channel" << '\n');
259 InteractionAvatar::restoreParticles();
260 return NULL;
261 }
262 }
263
264 /** Bias apply for this reaction in order to get the same
265 * ParticleBias for all stange particles.
266 * Can be reduced after because of the safeguard.
267 */
268 G4double bias_apply = 1.;
269 if(bias != 1.) bias_apply = Particle::getBiasFromVector(Particle::MergeVectorBias(particle1,particle2)) * bias;
270
271//// NN
272 if(particle1->isNucleon() && particle2->isNucleon()) {
273
274 G4double NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
275 G4double NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
276 G4double NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
277 G4double NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
278 G4double NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
279 G4double NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
280 G4double NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
281 G4double NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
282
283 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NNToNDelta(particle1, particle2) + CrossSections::NNToxPiNN(1,particle1, particle2)
284 + CrossSections::NNToxPiNN(2,particle1, particle2) + CrossSections::NNToxPiNN(3,particle1, particle2) + CrossSections::NNToxPiNN(4,particle1, particle2)
285 + CrossSections::NNToNNEtaExclu(particle1, particle2) + CrossSections::NNToNDeltaEta(particle1, particle2) + CrossSections::NNToNNEtaxPi(1,particle1, particle2)
286 + CrossSections::NNToNNEtaxPi(2,particle1, particle2) + CrossSections::NNToNNEtaxPi(3,particle1, particle2) + CrossSections::NNToNNEtaxPi(4,particle1, particle2)
287 + CrossSections::NNToNNOmegaExclu(particle1, particle2) + CrossSections::NNToNDeltaOmega(particle1, particle2) + CrossSections::NNToNNOmegaxPi(1,particle1, particle2)
288 + CrossSections::NNToNNOmegaxPi(2,particle1, particle2) + CrossSections::NNToNNOmegaxPi(3,particle1, particle2) + CrossSections::NNToNNOmegaxPi(4,particle1, particle2);
289 const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + NLKpiProductionCX + NSKpiProductionCX + NLK2piProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX)/bias_apply;
290
291 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
292
293 if(counterweight < 0.5) {
294 counterweight = 0.5;
295 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
296 NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
297 NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
298 NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
299 NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
300 NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
301 NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
302 NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
303 NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
304 }
305
306
307 const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
308 const G4double deltaProductionCX = CrossSections::NNToNDelta(particle1, particle2)*counterweight;
309 const G4double onePiProductionCX = CrossSections::NNToxPiNN(1,particle1, particle2)*counterweight;
310 const G4double twoPiProductionCX = CrossSections::NNToxPiNN(2,particle1, particle2)*counterweight;
311 const G4double threePiProductionCX = CrossSections::NNToxPiNN(3,particle1, particle2)*counterweight;
312 const G4double fourPiProductionCX = CrossSections::NNToxPiNN(4,particle1, particle2)*counterweight;
313
314 const G4double etaProductionCX = CrossSections::NNToNNEtaExclu(particle1, particle2)*counterweight;
315 const G4double etadeltaProductionCX = CrossSections::NNToNDeltaEta(particle1, particle2)*counterweight;
316 const G4double etaonePiProductionCX = CrossSections::NNToNNEtaxPi(1,particle1, particle2)*counterweight;
317 const G4double etatwoPiProductionCX = CrossSections::NNToNNEtaxPi(2,particle1, particle2)*counterweight;
318 const G4double etathreePiProductionCX = CrossSections::NNToNNEtaxPi(3,particle1, particle2)*counterweight;
319 const G4double etafourPiProductionCX = CrossSections::NNToNNEtaxPi(4,particle1, particle2)*counterweight;
320 const G4double omegaProductionCX = CrossSections::NNToNNOmegaExclu(particle1, particle2)*counterweight;
321 const G4double omegadeltaProductionCX = CrossSections::NNToNDeltaOmega(particle1, particle2)*counterweight;
322 const G4double omegaonePiProductionCX = CrossSections::NNToNNOmegaxPi(1,particle1, particle2)*counterweight;
323 const G4double omegatwoPiProductionCX = CrossSections::NNToNNOmegaxPi(2,particle1, particle2)*counterweight;
324 const G4double omegathreePiProductionCX = CrossSections::NNToNNOmegaxPi(3,particle1, particle2)*counterweight;
325 const G4double omegafourPiProductionCX = CrossSections::NNToNNOmegaxPi(4,particle1, particle2)*counterweight;
326
327 const G4double totCX=CrossSections::total(particle1, particle2);
328
329// assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-fourPiProductionCX-NLKProductionCX-NSKProductionCX-NLKpiProductionCX-NSKpiProductionCX-NLK2piProductionCX-NSK2piProductionCX-NNKKbProductionCX-NNMissingCX-etaProductionCX-etadeltaProductionCX-etaonePiProductionCX-etatwoPiProductionCX-etathreePiProductionCX-etafourPiProductionCX-omegaProductionCX-omegadeltaProductionCX-omegaonePiProductionCX-omegatwoPiProductionCX-omegathreePiProductionCX-omegafourPiProductionCX) < 0.5);
330
331 const G4double rChannel=Random::shoot() * totCX;
332
333 if(elasticCX > rChannel) {
334// Elastic NN channel
335 isElastic = true;
336 weight = counterweight;
337 if (theNucleus->getStore()->getBook().getAcceptedCollisions() == 0 &&
338 theNucleus->getStore()->getBook().getAcceptedSrcCollisions() == 0 &&
339 ParticleTable::getsrcPairConfig()) {
340 INCL_DEBUG("NN-SRC interaction: elastic channel chosen" << '\n');
341 ParticleList &inside = theNucleus->getStore()->getParticlesforSrc();
342 G4int zz = theNucleus->getZ();
343 generateSrcPairsMethod(inside, theNucleus->getA() - zz, zz, particle1,
344 particle2);
345 if ((particle1->getSrcPair() > 0 || particle2->getSrcPair() > 0)) {
346 return new ElasticChannel(particle1, particle2, theNucleus);
347 } else {
348 INCL_DEBUG("NN interaction: elastic channel chosen" << '\n');
349 theNucleus->resetSrc();
350 return new ElasticChannel(particle1, particle2);
351 }
352 } else {
353 INCL_DEBUG("NN interaction: elastic channel chosen" << '\n');
354 return new ElasticChannel(particle1, particle2);
355 }
356 } else if((elasticCX + deltaProductionCX) > rChannel) {
357 isElastic = false;
358// NN -> N Delta channel is chosen
359 weight = counterweight;
360 if (theNucleus->getStore()->getBook().getAcceptedCollisions() == 0 &&
361 theNucleus->getStore()->getBook().getAcceptedSrcCollisions() == 0 &&
362 ParticleTable::getsrcPairConfig()) {
363 INCL_DEBUG("NN-SRC interaction: Delta channel chosen" << '\n');
364 ParticleList &inside = theNucleus->getStore()->getParticlesforSrc();
365 G4int zz = theNucleus->getZ();
366 generateSrcPairsMethod(inside, theNucleus->getA() - zz, zz, particle1,
367 particle2);
368
369 if ((particle1->getSrcPair() > 0 || particle2->getSrcPair() > 0)) {
370 return new DeltaProductionChannel(particle1, particle2, theNucleus);
371 } else {
372 INCL_DEBUG("NN interaction: Delta channel chosen" << '\n');
373 theNucleus->resetSrc();
374 return new DeltaProductionChannel(particle1, particle2);
375 }
376 } else {
377 INCL_DEBUG("NN interaction: Delta channel chosen" << '\n');
378 return new DeltaProductionChannel(particle1, particle2);
379 }
380 } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
381 isElastic = false;
382// NN -> PiNN channel is chosen
383 INCL_DEBUG("NN interaction: one Pion channel chosen" << '\n');
384 weight = counterweight;
385 return new NNToMultiPionsChannel(1,particle1, particle2);
386 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
387 isElastic = false;
388// NN -> 2PiNN channel is chosen
389 INCL_DEBUG("NN interaction: two Pions channel chosen" << '\n');
390 weight = counterweight;
391 return new NNToMultiPionsChannel(2,particle1, particle2);
392 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
393 isElastic = false;
394// NN -> 3PiNN channel is chosen
395 INCL_DEBUG("NN interaction: three Pions channel chosen" << '\n');
396 weight = counterweight;
397 return new NNToMultiPionsChannel(3,particle1, particle2);
398 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX > rChannel) {
399 isElastic = false;
400// NN -> 4PiNN channel is chosen
401 INCL_DEBUG("NN interaction: four Pions channel chosen" << '\n');
402 weight = counterweight;
403 return new NNToMultiPionsChannel(4,particle1, particle2);
404 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
405 + etaProductionCX > rChannel) {
406 isElastic = false;
407// NN -> NNEta channel is chosen
408 INCL_DEBUG("NN interaction: Eta channel chosen" << '\n');
409 weight = counterweight;
410 return new NNToNNEtaChannel(particle1, particle2);
411 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
412 + etaProductionCX + etadeltaProductionCX > rChannel) {
413 isElastic = false;
414// NN -> N Delta Eta channel is chosen
415 INCL_DEBUG("NN interaction: Delta Eta channel chosen" << '\n');
416 weight = counterweight;
417 return new NDeltaEtaProductionChannel(particle1, particle2);
418 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
419 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX > rChannel) {
420 isElastic = false;
421// NN -> EtaPiNN channel is chosen
422 INCL_DEBUG("NN interaction: Eta + one Pion channel chosen" << '\n');
423 weight = counterweight;
424 return new NNEtaToMultiPionsChannel(1,particle1, particle2);
425 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
426 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX > rChannel) {
427 isElastic = false;
428// NN -> Eta2PiNN channel is chosen
429 INCL_DEBUG("NN interaction: Eta + two Pions channel chosen" << '\n');
430 weight = counterweight;
431 return new NNEtaToMultiPionsChannel(2,particle1, particle2);
432 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
433 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX > rChannel) {
434 isElastic = false;
435// NN -> Eta3PiNN channel is chosen
436 INCL_DEBUG("NN interaction: Eta + three Pions channel chosen" << '\n');
437 weight = counterweight;
438 return new NNEtaToMultiPionsChannel(3,particle1, particle2);
439 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
440 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX > rChannel) {
441 isElastic = false;
442// NN -> Eta4PiNN channel is chosen
443 INCL_DEBUG("NN interaction: Eta + four Pions channel chosen" << '\n');
444 weight = counterweight;
445 return new NNEtaToMultiPionsChannel(4,particle1, particle2);
446 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
447 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
448 + omegaProductionCX > rChannel) {
449 isElastic = false;
450// NN -> NNOmega channel is chosen
451 INCL_DEBUG("NN interaction: Omega channel chosen" << '\n');
452 weight = counterweight;
453 return new NNToNNOmegaChannel(particle1, particle2);
454 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
455 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
456 + omegaProductionCX + omegadeltaProductionCX > rChannel) {
457 isElastic = false;
458// NN -> N Delta Omega channel is chosen
459 INCL_DEBUG("NN interaction: Delta Omega channel chosen" << '\n');
460 weight = counterweight;
461 return new NDeltaOmegaProductionChannel(particle1, particle2);
462 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
463 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
464 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX > rChannel) {
465 isElastic = false;
466// NN -> OmegaPiNN channel is chosen
467 INCL_DEBUG("NN interaction: Omega + one Pion channel chosen" << '\n');
468 weight = counterweight;
469 return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
470 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
471 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
472 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX > rChannel) {
473 isElastic = false;
474// NN -> Omega2PiNN channel is chosen
475 INCL_DEBUG("NN interaction: Omega + two Pions channel chosen" << '\n');
476 weight = counterweight;
477 return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
478 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
479 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
480 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX > rChannel) {
481 isElastic = false;
482// NN -> Omega3PiNN channel is chosen
483 INCL_DEBUG("NN interaction: Omega + three Pions channel chosen" << '\n');
484 weight = counterweight;
485 return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
486 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
487 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
488 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX > rChannel) {
489 isElastic = false;
490// NN -> Omega4PiNN channel is chosen
491 INCL_DEBUG("NN interaction: Omega + four Pions channel chosen" << '\n');
492 weight = counterweight;
493 return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
494 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
495 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
496 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
497 + NLKProductionCX > rChannel) {
498 isElastic = false;
499 isStrangeProduction = true;
500// NN -> NLK channel is chosen
501 INCL_DEBUG("NN interaction: NLK channel chosen" << '\n');
502 weight = bias_apply;
503 return new NNToNLKChannel(particle1, particle2);
504 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
505 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
506 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
507 + NLKProductionCX + NLKpiProductionCX > rChannel) {
508 isElastic = false;
509 isStrangeProduction = true;
510// NN -> NLKpi channel is chosen
511 INCL_DEBUG("NN interaction: NLKpi channel chosen" << '\n');
512 weight = bias_apply;
513 return new NNToNLKpiChannel(particle1, particle2);
514 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
515 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
516 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
517 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX > rChannel) {
518 isElastic = false;
519 isStrangeProduction = true;
520// NN -> NLK2pi channel is chosen
521 INCL_DEBUG("NN interaction: NLK2pi channel chosen" << '\n');
522 weight = bias_apply;
523 return new NNToNLK2piChannel(particle1, particle2);
524 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
525 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
526 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
527 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX > rChannel) {
528 isElastic = false;
529 isStrangeProduction = true;
530// NN -> NSK channel is chosen
531 INCL_DEBUG("NN interaction: NSK channel chosen" << '\n');
532 weight = bias_apply;
533 return new NNToNSKChannel(particle1, particle2);
534 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
535 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
536 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
537 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX > rChannel) {
538 isElastic = false;
539 isStrangeProduction = true;
540// NN -> NSKpi channel is chosen
541 INCL_DEBUG("NN interaction: NSKpi channel chosen" << '\n');
542 weight = bias_apply;
543 return new NNToNSKpiChannel(particle1, particle2);
544 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
545 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
546 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
547 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX > rChannel) {
548 isElastic = false;
549 isStrangeProduction = true;
550// NN -> NSK2pi channel is chosen
551 INCL_DEBUG("NN interaction: NSK2pi channel chosen" << '\n');
552 weight = bias_apply;
553 return new NNToNSK2piChannel(particle1, particle2);
554 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
555 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
556 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
557 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX > rChannel) {
558 isElastic = false;
559 isStrangeProduction = true;
560// NN -> NNKKb channel is chosen
561 INCL_DEBUG("NN interaction: NNKKb channel chosen" << '\n');
562 weight = bias_apply;
563 return new NNToNNKKbChannel(particle1, particle2);
564 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
565 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
566 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
567 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX> rChannel) {
568 isElastic = false;
569 isStrangeProduction = true;
570// NN -> Missing Strangeness channel is chosen
571 INCL_DEBUG("NN interaction: Missing Strangeness channel chosen" << '\n');
572 weight = bias_apply;
573 return new NNToMissingStrangenessChannel(particle1, particle2);
574 } else {
575 INCL_WARN("inconsistency within the NN Cross Sections (sum!=inelastic)" << '\n');
576 if(NNMissingCX>0.) {
577 INCL_WARN("Returning an Missing Strangeness channel" << '\n');
578 weight = bias_apply;
579 isElastic = false;
580 isStrangeProduction = true;
581 return new NNToNNKKbChannel(particle1, particle2);
582 } else if(NNKKbProductionCX>0.) {
583 INCL_WARN("Returning an NNKKb channel" << '\n');
584 weight = bias_apply;
585 isElastic = false;
586 isStrangeProduction = true;
587 return new NNToNNKKbChannel(particle1, particle2);
588 } else if(NSK2piProductionCX>0.) {
589 INCL_WARN("Returning an NSK2pi channel" << '\n');
590 weight = bias_apply;
591 isElastic = false;
592 isStrangeProduction = true;
593 return new NNToNSK2piChannel(particle1, particle2);
594 } else if(NSKpiProductionCX>0.) {
595 INCL_WARN("Returning an NSKpi channel" << '\n');
596 weight = bias_apply;
597 isElastic = false;
598 isStrangeProduction = true;
599 return new NNToNSKpiChannel(particle1, particle2);
600 } else if(NSKProductionCX>0.) {
601 INCL_WARN("Returning an NSK channel" << '\n');
602 weight = bias_apply;
603 isElastic = false;
604 isStrangeProduction = true;
605 return new NNToNSKChannel(particle1, particle2);
606 } else if(NLK2piProductionCX>0.) {
607 INCL_WARN("Returning an NLK2pi channel" << '\n');
608 weight = bias_apply;
609 isElastic = false;
610 isStrangeProduction = true;
611 return new NNToNLK2piChannel(particle1, particle2);
612 } else if(NLKpiProductionCX>0.) {
613 INCL_WARN("Returning an NLKpi channel" << '\n');
614 weight = bias_apply;
615 isElastic = false;
616 isStrangeProduction = true;
617 return new NNToNLKpiChannel(particle1, particle2);
618 } else if(NLKProductionCX>0.) {
619 INCL_WARN("Returning an NLK channel" << '\n');
620 weight = bias_apply;
621 isElastic = false;
622 isStrangeProduction = true;
623 return new NNToNLKChannel(particle1, particle2);
624 } else if(omegafourPiProductionCX>0.) {
625 INCL_WARN("Returning an Omega + four Pions channel" << '\n');
626 weight = counterweight;
627 isElastic = false;
628 return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
629 } else if(omegathreePiProductionCX>0.) {
630 INCL_WARN("Returning an Omega + three Pions channel" << '\n');
631 weight = counterweight;
632 isElastic = false;
633 return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
634 } else if(omegatwoPiProductionCX>0.) {
635 INCL_WARN("Returning an Omega + two Pions channel" << '\n');
636 weight = counterweight;
637 isElastic = false;
638 return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
639 } else if(omegaonePiProductionCX>0.) {
640 INCL_WARN("Returning an Omega + one Pion channel" << '\n');
641 weight = counterweight;
642 isElastic = false;
643 return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
644 } else if(omegadeltaProductionCX>0.) {
645 INCL_WARN("Returning an Omega + Delta channel" << '\n');
646 weight = counterweight;
647 isElastic = false;
648 return new NDeltaOmegaProductionChannel(particle1, particle2);
649 } else if(omegaProductionCX>0.) {
650 INCL_WARN("Returning an Omega channel" << '\n');
651 weight = counterweight;
652 isElastic = false;
653 return new NNToNNOmegaChannel(particle1, particle2);
654 } else if(etafourPiProductionCX>0.) {
655 INCL_WARN("Returning an Eta + four Pions channel" << '\n');
656 weight = counterweight;
657 isElastic = false;
658 return new NNEtaToMultiPionsChannel(4,particle1, particle2);
659 } else if(etathreePiProductionCX>0.) {
660 INCL_WARN("Returning an Eta + threev channel" << '\n');
661 weight = counterweight;
662 isElastic = false;
663 return new NNEtaToMultiPionsChannel(3,particle1, particle2);
664 } else if(etatwoPiProductionCX>0.) {
665 INCL_WARN("Returning an Eta + two Pions channel" << '\n');
666 weight = counterweight;
667 isElastic = false;
668 return new NNEtaToMultiPionsChannel(2,particle1, particle2);
669 } else if(etaonePiProductionCX>0.) {
670 INCL_WARN("Returning an Eta + one Pion channel" << '\n');
671 weight = counterweight;
672 isElastic = false;
673 return new NNEtaToMultiPionsChannel(1,particle1, particle2);
674 } else if(etadeltaProductionCX>0.) {
675 INCL_WARN("Returning an Eta + Delta channel" << '\n');
676 weight = counterweight;
677 isElastic = false;
678 return new NDeltaEtaProductionChannel(particle1, particle2);
679 } else if(etaProductionCX>0.) {
680 INCL_WARN("Returning an Eta channel" << '\n');
681 weight = counterweight;
682 isElastic = false;
683 return new NNToNNEtaChannel(particle1, particle2);
684 } else if(fourPiProductionCX>0.) {
685 INCL_WARN("Returning a 4pi channel" << '\n');
686 weight = counterweight;
687 isElastic = false;
688 return new NNToMultiPionsChannel(4,particle1, particle2);
689 } else if(threePiProductionCX>0.) {
690 INCL_WARN("Returning a 3pi channel" << '\n');
691 weight = counterweight;
692 isElastic = false;
693 return new NNToMultiPionsChannel(3,particle1, particle2);
694 } else if(twoPiProductionCX>0.) {
695 INCL_WARN("Returning a 2pi channel" << '\n');
696 weight = counterweight;
697 isElastic = false;
698 return new NNToMultiPionsChannel(2,particle1, particle2);
699 } else if(onePiProductionCX>0.) {
700 INCL_WARN("Returning a 1pi channel" << '\n');
701 weight = counterweight;
702 isElastic = false;
703 return new NNToMultiPionsChannel(1,particle1, particle2);
704 } else if(deltaProductionCX>0.) {
705 INCL_WARN("Returning a delta-production channel" << '\n');
706 weight = counterweight;
707 isElastic = false;
708 return new DeltaProductionChannel(particle1, particle2);
709 } else {
710 INCL_WARN("Returning an elastic channel" << '\n');
711 weight = counterweight;
712 isElastic = true;
713 return new ElasticChannel(particle1, particle2);
714 }
715 }
716
717//// NDelta
718 }
719 else if((particle1->isNucleon() && particle2->isDelta()) ||
720 (particle1->isDelta() && particle2->isNucleon())) {
721
722 G4double NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
723 G4double NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
724 G4double DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
725 G4double DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
726 G4double NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
727
728 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NDeltaToNN(particle1, particle2);
729 const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX)/bias_apply;
730
731 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
732
733 if(counterweight < 0.5){
734 counterweight = 0.5;
735 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
736
737 NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
738 NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
739 DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
740 DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
741 NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
742 }
743
744 G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
745 G4double recombinationCX = CrossSections::NDeltaToNN(particle1, particle2)*counterweight;
746
747 const G4double rChannel=Random::shoot() * (StrangenessProdCX + UnStrangeProdCX);
748
749 if(elasticCX > rChannel) {
750 isElastic = true;
751// Elastic N Delta channel
752 INCL_DEBUG("NDelta interaction: elastic channel chosen" << '\n');
753 weight = counterweight;
754 return new ElasticChannel(particle1, particle2);
755 } else if (elasticCX + recombinationCX > rChannel){
756 isElastic = false;
757// Recombination
758// NDelta -> NN channel is chosen
759 INCL_DEBUG("NDelta interaction: recombination channel chosen" << '\n');
760 weight = counterweight;
761 return new RecombinationChannel(particle1, particle2);
762 } else if (elasticCX + recombinationCX + NLKProductionCX > rChannel){
763 isElastic = false;
764 isStrangeProduction = true;
765// NDelta -> NLK channel is chosen
766 INCL_DEBUG("NDelta interaction: NLK channel chosen" << '\n');
767 weight = bias_apply;
768 return new NDeltaToNLKChannel(particle1, particle2);
769 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX > rChannel){
770 isElastic = false;
771 isStrangeProduction = true;
772// NDelta -> NSK channel is chosen
773 INCL_DEBUG("NDelta interaction: NSK channel chosen" << '\n');
774 weight = bias_apply;
775 return new NDeltaToNSKChannel(particle1, particle2);
776 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX > rChannel){
777 isElastic = false;
778 isStrangeProduction = true;
779// NDelta -> DeltaLK channel is chosen
780 INCL_DEBUG("NDelta interaction: DeltaLK channel chosen" << '\n');
781 weight = bias_apply;
782 return new NDeltaToDeltaLKChannel(particle1, particle2);
783 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX > rChannel){
784 isElastic = false;
785 isStrangeProduction = true;
786// NDelta -> DeltaSK channel is chosen
787 INCL_DEBUG("NDelta interaction: DeltaSK channel chosen" << '\n');
788 weight = bias_apply;
789 return new NDeltaToDeltaSKChannel(particle1, particle2);
790 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX > rChannel){
791 isElastic = false;
792 isStrangeProduction = true;
793// NDelta -> NNKKb channel is chosen
794 INCL_DEBUG("NDelta interaction: NNKKb channel chosen" << '\n');
795 weight = bias_apply;
796 return new NDeltaToNNKKbChannel(particle1, particle2);
797 }
798 else{
799 INCL_ERROR("rChannel > (StrangenessProdCX + UnStrangeProdCX) in NDelta interaction: return an elastic channel" << '\n');
800 weight = counterweight;
801 isElastic = true;
802 return new ElasticChannel(particle1, particle2);
803 }
804
805//// DeltaDelta
806 } else if(particle1->isDelta() && particle2->isDelta()) {
807 isElastic = true;
808 INCL_DEBUG("DeltaDelta interaction: elastic channel chosen" << '\n');
809 return new ElasticChannel(particle1, particle2);
810
811//// PiN
812 } else if(isPiN) {
813
814 G4double LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
815 G4double SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
816 G4double LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
817 G4double SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
818 G4double LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
819 G4double SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
820 G4double NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
821 G4double MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
822
823 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::piNToDelta(particle1, particle2)
824 + CrossSections::piNToxPiN(2,particle1, particle2) + CrossSections::piNToxPiN(3,particle1, particle2) + CrossSections::piNToxPiN(4,particle1, particle2)
825 + CrossSections::piNToEtaN(particle1, particle2) + CrossSections::piNToOmegaN(particle1, particle2);
826 const G4double StrangenessProdCX = (LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX)/bias_apply;
827
828 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
829
830 if(counterweight < 0.5) {
831 counterweight = 0.5;
832 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
833 LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
834 SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
835 LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
836 SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
837 LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
838 SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
839 NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
840 MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
841 }
842
843
844 const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
845 const G4double deltaProductionCX = CrossSections::piNToDelta(particle1, particle2)*counterweight;
846 const G4double onePiProductionCX = CrossSections::piNToxPiN(2,particle1, particle2)*counterweight;
847 const G4double twoPiProductionCX = CrossSections::piNToxPiN(3,particle1, particle2)*counterweight;
848 const G4double threePiProductionCX = CrossSections::piNToxPiN(4,particle1, particle2)*counterweight;
849 const G4double etaProductionCX = CrossSections::piNToEtaN(particle1, particle2)*counterweight;
850 const G4double omegaProductionCX = CrossSections::piNToOmegaN(particle1, particle2)*counterweight;
851
852 const G4double totCX=CrossSections::total(particle1, particle2);
853
854// assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-etaProductionCX-omegaProductionCX-LKProdCX-SKProdCX-LKpiProdCX-SKpiProdCX-LK2piProdCX-SK2piProdCX-NKKbProdCX-MissingCX) < 0.15);
855
856 const G4double rChannel=Random::shoot() * totCX;
857
858 if(elasticCX > rChannel) {
859 isElastic = true;
860// Elastic PiN channel
861 INCL_DEBUG("PiN interaction: elastic channel chosen" << '\n');
862 weight = counterweight;
863 return new PiNElasticChannel(particle1, particle2);
864 } else if(elasticCX + deltaProductionCX > rChannel) {
865 isElastic = false;
866// PiN -> Delta channel is chosen
867 INCL_DEBUG("PiN interaction: Delta channel chosen" << '\n');
868 weight = counterweight;
869 return new PiNToDeltaChannel(particle1, particle2);
870 } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
871 isElastic = false;
872// PiN -> PiNPi channel is chosen
873 INCL_DEBUG("PiN interaction: one Pion channel chosen" << '\n');
874 weight = counterweight;
875 return new PiNToMultiPionsChannel(2,particle1, particle2);
876 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
877 isElastic = false;
878// PiN -> PiN2Pi channel is chosen
879 INCL_DEBUG("PiN interaction: two Pions channel chosen" << '\n');
880 weight = counterweight;
881 return new PiNToMultiPionsChannel(3,particle1, particle2);
882 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
883 isElastic = false;
884// PiN -> PiN3Pi channel is chosen
885 INCL_DEBUG("PiN interaction: three Pions channel chosen" << '\n');
886 weight = counterweight;
887 return new PiNToMultiPionsChannel(4,particle1, particle2);
888 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX > rChannel) {
889 isElastic = false;
890// PiN -> EtaN channel is chosen
891 INCL_DEBUG("PiN interaction: Eta channel chosen" << '\n');
892 weight = counterweight;
893 return new PiNToEtaChannel(particle1, particle2);
894 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX > rChannel) {
895 isElastic = false;
896// PiN -> OmegaN channel is chosen
897 INCL_DEBUG("PiN interaction: Omega channel chosen" << '\n');
898 weight = counterweight;
899 return new PiNToOmegaChannel(particle1, particle2);
900 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
901 + LKProdCX > rChannel) {
902 isElastic = false;
903 isStrangeProduction = true;
904// PiN -> LK channel is chosen
905 INCL_DEBUG("PiN interaction: LK channel chosen" << '\n');
906 weight = bias_apply;
907 return new NpiToLKChannel(particle1, particle2);
908 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
909 + LKProdCX + SKProdCX > rChannel) {
910 isElastic = false;
911 isStrangeProduction = true;
912// PiN -> SK channel is chosen
913 INCL_DEBUG("PiN interaction: SK channel chosen" << '\n');
914 weight = bias_apply;
915 return new NpiToSKChannel(particle1, particle2);
916 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
917 + LKProdCX + SKProdCX + LKpiProdCX > rChannel) {
918 isElastic = false;
919 isStrangeProduction = true;
920// PiN -> LKpi channel is chosen
921 INCL_DEBUG("PiN interaction: LKpi channel chosen" << '\n');
922 weight = bias_apply;
923 return new NpiToLKpiChannel(particle1, particle2);
924 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
925 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX > rChannel) {
926 isElastic = false;
927 isStrangeProduction = true;
928// PiN -> SKpi channel is chosen
929 INCL_DEBUG("PiN interaction: SKpi channel chosen" << '\n');
930 weight = bias_apply;
931 return new NpiToSKpiChannel(particle1, particle2);
932 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
933 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX > rChannel) {
934 isElastic = false;
935 isStrangeProduction = true;
936// PiN -> LK2pi channel is chosen
937 INCL_DEBUG("PiN interaction: LK2pi channel chosen" << '\n');
938 weight = bias_apply;
939 return new NpiToLK2piChannel(particle1, particle2);
940 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
941 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX > rChannel) {
942 isElastic = false;
943 isStrangeProduction = true;
944// PiN -> SK2pi channel is chosen
945 INCL_DEBUG("PiN interaction: SK2pi channel chosen" << '\n');
946 weight = bias_apply;
947 return new NpiToSK2piChannel(particle1, particle2);
948 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
949 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX > rChannel) {
950 isElastic = false;
951 isStrangeProduction = true;
952// PiN -> NKKb channel is chosen
953 INCL_DEBUG("PiN interaction: NKKb channel chosen" << '\n');
954 weight = bias_apply;
955 return new NpiToNKKbChannel(particle1, particle2);
956 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
957 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX> rChannel) {
958 isElastic = false;
959 isStrangeProduction = true;
960// PiN -> Missinge Strangeness channel is chosen
961 INCL_DEBUG("PiN interaction: Missinge Strangeness channel chosen" << '\n');
962 weight = bias_apply;
963 return new NpiToMissingStrangenessChannel(particle1, particle2);
964 }
965 else {
966 INCL_WARN("inconsistency within the PiN Cross Sections (sum!=inelastic)" << '\n');
967 if(MissingCX>0.) {
968 INCL_WARN("Returning a Missinge Strangeness channel" << '\n');
969 weight = bias_apply;
970 isElastic = false;
971 isStrangeProduction = true;
972 return new NpiToMissingStrangenessChannel(particle1, particle2);
973 } else if(NKKbProdCX>0.) {
974 INCL_WARN("Returning a NKKb channel" << '\n');
975 weight = bias_apply;
976 isElastic = false;
977 isStrangeProduction = true;
978 return new NpiToNKKbChannel(particle1, particle2);
979 } else if(SK2piProdCX>0.) {
980 INCL_WARN("Returning a SK2pi channel" << '\n');
981 weight = bias_apply;
982 isElastic = false;
983 isStrangeProduction = true;
984 return new NpiToSK2piChannel(particle1, particle2);
985 } else if(LK2piProdCX>0.) {
986 INCL_WARN("Returning a LK2pi channel" << '\n');
987 weight = bias_apply;
988 isElastic = false;
989 isStrangeProduction = true;
990 return new NpiToLK2piChannel(particle1, particle2);
991 } else if(SKpiProdCX>0.) {
992 INCL_WARN("Returning a SKpi channel" << '\n');
993 weight = bias_apply;
994 isElastic = false;
995 isStrangeProduction = true;
996 return new NpiToSKpiChannel(particle1, particle2);
997 } else if(LKpiProdCX>0.) {
998 INCL_WARN("Returning a LKpi channel" << '\n');
999 weight = bias_apply;
1000 isElastic = false;
1001 isStrangeProduction = true;
1002 return new NpiToLKpiChannel(particle1, particle2);
1003 } else if(SKProdCX>0.) {
1004 INCL_WARN("Returning a SK channel" << '\n');
1005 weight = bias_apply;
1006 isElastic = false;
1007 isStrangeProduction = true;
1008 return new NpiToSKChannel(particle1, particle2);
1009 } else if(LKProdCX>0.) {
1010 INCL_WARN("Returning a LK channel" << '\n');
1011 weight = bias_apply;
1012 isElastic = false;
1013 isStrangeProduction = true;
1014 return new NpiToLKChannel(particle1, particle2);
1015 } else if(omegaProductionCX>0.) {
1016 INCL_WARN("Returning a Omega channel" << '\n');
1017 weight = counterweight;
1018 isElastic = false;
1019 return new PiNToOmegaChannel(particle1, particle2);
1020 } else if(etaProductionCX>0.) {
1021 INCL_WARN("Returning a Eta channel" << '\n');
1022 weight = counterweight;
1023 isElastic = false;
1024 return new PiNToEtaChannel(particle1, particle2);
1025 } else if(threePiProductionCX>0.) {
1026 INCL_WARN("Returning a 3pi channel" << '\n');
1027 weight = counterweight;
1028 isElastic = false;
1029 return new PiNToMultiPionsChannel(4,particle1, particle2);
1030 } else if(twoPiProductionCX>0.) {
1031 INCL_WARN("Returning a 2pi channel" << '\n');
1032 weight = counterweight;
1033 isElastic = false;
1034 return new PiNToMultiPionsChannel(3,particle1, particle2);
1035 } else if(onePiProductionCX>0.) {
1036 INCL_WARN("Returning a 1pi channel" << '\n');
1037 weight = counterweight;
1038 isElastic = false;
1039 return new PiNToMultiPionsChannel(2,particle1, particle2);
1040 } else if(deltaProductionCX>0.) {
1041 INCL_WARN("Returning a delta-production channel" << '\n');
1042 weight = counterweight;
1043 isElastic = false;
1044 return new PiNToDeltaChannel(particle1, particle2);
1045 } else {
1046 INCL_WARN("Returning an elastic channel" << '\n');
1047 weight = counterweight;
1048 isElastic = true;
1049 return new PiNElasticChannel(particle1, particle2);
1050 }
1051 }
1052 } else if ((particle1->isNucleon() && particle2->isEta()) || (particle2->isNucleon() && particle1->isEta())) {
1053//// EtaN
1054
1055 const G4double elasticCX = CrossSections::elastic(particle1, particle2);
1056 const G4double onePiProductionCX = CrossSections::etaNToPiN(particle1, particle2);
1057 const G4double twoPiProductionCX = CrossSections::etaNToPiPiN(particle1, particle2);
1058 const G4double LKProductionCX = CrossSections::etaNToLK(particle1, particle2);
1059 const G4double SKProductionCX = CrossSections::etaNToSK(particle1, particle2);
1060 const G4double totCX=CrossSections::total(particle1, particle2);
1061// assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)-LKProductionCX-SKProductionCX<1.);
1062
1063 const G4double rChannel=Random::shoot() * totCX;
1064
1065 if(elasticCX > rChannel) {
1066// Elastic EtaN channel
1067 isElastic = true;
1068 INCL_DEBUG("EtaN interaction: elastic channel chosen" << '\n');
1069 return new EtaNElasticChannel(particle1, particle2);
1070 } else if(elasticCX + onePiProductionCX > rChannel) {
1071 isElastic = false;
1072// EtaN -> EtaPiN channel is chosen
1073 INCL_DEBUG("EtaN interaction: PiN channel chosen" << '\n');
1074 return new EtaNToPiNChannel(particle1, particle2);
1075 } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
1076 isElastic = false;
1077// EtaN -> EtaPiPiN channel is chosen
1078 INCL_DEBUG("EtaN interaction: PiPiN channel chosen" << '\n');
1079 return new EtaNToPiPiNChannel(particle1, particle2);
1080 } else if(elasticCX + onePiProductionCX + twoPiProductionCX + LKProductionCX > rChannel) {
1081 isElastic = false;
1082// EtaN -> LK channel is chosen
1083 INCL_DEBUG("EtaN interaction: LK channel chosen" << '\n');
1084 return new EtaOrOmegaNToLKChannel(particle1, particle2);
1085 } else if(elasticCX + onePiProductionCX + twoPiProductionCX + LKProductionCX + SKProductionCX > rChannel) {
1086 isElastic = false;
1087// EtaN -> SK channel is chosen
1088 INCL_DEBUG("EtaN interaction: SK channel chosen" << '\n');
1089 return new EtaOrOmegaNToSKChannel(particle1, particle2);
1090 }
1091
1092 else {
1093 INCL_WARN("inconsistency within the EtaN Cross Sections (sum!=inelastic)" << '\n');
1094 if(SKProductionCX>0.) {
1095 INCL_WARN("Returning a SK channel" << '\n');
1096 isElastic = false;
1097 return new EtaOrOmegaNToSKChannel(particle1, particle2);
1098 } else if(LKProductionCX>0.) {
1099 INCL_WARN("Returning a LK channel" << '\n');
1100 isElastic = false;
1101 return new EtaOrOmegaNToLKChannel(particle1, particle2);
1102 } else if(twoPiProductionCX>0.) {
1103 INCL_WARN("Returning a PiPiN channel" << '\n');
1104 isElastic = false;
1105 return new EtaNToPiPiNChannel(particle1, particle2);
1106 } else if(onePiProductionCX>0.) {
1107 INCL_WARN("Returning a PiN channel" << '\n');
1108 isElastic = false;
1109 return new EtaNToPiNChannel(particle1, particle2);
1110 } else {
1111 INCL_WARN("Returning an elastic channel" << '\n');
1112 isElastic = true;
1113 return new EtaNElasticChannel(particle1, particle2);
1114 }
1115 }
1116
1117 } else if ((particle1->isNucleon() && particle2->isOmega()) || (particle2->isNucleon() && particle1->isOmega())) {
1118//// OmegaN
1119
1120 const G4double elasticCX = CrossSections::elastic(particle1, particle2);
1121 const G4double onePiProductionCX = CrossSections::omegaNToPiN(particle1, particle2);
1122 const G4double twoPiProductionCX = CrossSections::omegaNToPiPiN(particle1, particle2);
1123 const G4double LKProductionCX = CrossSections::omegaNToLK(particle1, particle2);
1124 const G4double SKProductionCX = CrossSections::omegaNToSK(particle1, particle2);
1125 const G4double totCX=CrossSections::total(particle1, particle2);
1126// assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX + LKProductionCX + SKProductionCX)<1.);
1127
1128 const G4double rChannel=Random::shoot() * totCX;
1129
1130 if(elasticCX > rChannel) {
1131// Elastic OmegaN channel
1132 isElastic = true;
1133 INCL_DEBUG("OmegaN interaction: elastic channel chosen" << '\n');
1134 return new OmegaNElasticChannel(particle1, particle2);
1135 } else if(elasticCX + onePiProductionCX > rChannel) {
1136 isElastic = false;
1137// OmegaN -> PiN channel is chosen
1138 INCL_DEBUG("OmegaN interaction: PiN channel chosen" << '\n');
1139 return new OmegaNToPiNChannel(particle1, particle2);
1140 } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
1141 isElastic = false;
1142// OmegaN -> PiPiN channel is chosen
1143 INCL_DEBUG("OmegaN interaction: PiPiN channel chosen" << '\n');
1144 return new OmegaNToPiPiNChannel(particle1, particle2);
1145 } else if(elasticCX + onePiProductionCX + twoPiProductionCX + LKProductionCX > rChannel) {
1146 isElastic = false;
1147// OmegaN -> LK channel is chosen
1148 INCL_DEBUG("EtaN interaction: LK channel chosen" << '\n');
1149 return new EtaOrOmegaNToLKChannel(particle1, particle2);
1150 } else if(elasticCX + onePiProductionCX + twoPiProductionCX + LKProductionCX + SKProductionCX > rChannel) {
1151 isElastic = false;
1152// OmegaN -> SK channel is chosen
1153 INCL_DEBUG("EtaN interaction: SK channel chosen" << '\n');
1154 return new EtaOrOmegaNToSKChannel(particle1, particle2);
1155 }
1156 else {
1157 INCL_WARN("inconsistency within the OmegaN Cross Sections (sum!=inelastic)" << '\n');
1158 if(SKProductionCX>0.) {
1159 INCL_WARN("Returning a SK channel" << '\n');
1160 isElastic = false;
1161 return new EtaOrOmegaNToSKChannel(particle1, particle2);
1162 } else if(LKProductionCX>0.) {
1163 INCL_WARN("Returning a LK channel" << '\n');
1164 isElastic = false;
1165 return new EtaOrOmegaNToLKChannel(particle1, particle2);
1166 } else if(twoPiProductionCX>0.) {
1167 INCL_WARN("Returning a PiPiN channel" << '\n');
1168 isElastic = false;
1169 return new OmegaNToPiPiNChannel(particle1, particle2);
1170 } else if(onePiProductionCX>0.) {
1171 INCL_WARN("Returning a PiN channel" << '\n');
1172 isElastic = false;
1173 return new OmegaNToPiNChannel(particle1, particle2);
1174 } else {
1175 INCL_WARN("Returning an elastic channel" << '\n');
1176 isElastic = true;
1177 return new OmegaNElasticChannel(particle1, particle2);
1178 }
1179 }
1180 } else if ((particle1->isNucleon() && particle2->isKaon()) || (particle2->isNucleon() && particle1->isKaon())) {
1181//// KN
1182 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1183 const G4double quasielasticCX = CrossSections::NKToNK(particle1,particle2);
1184 const G4double NKToNKpiCX = CrossSections::NKToNKpi(particle1,particle2);
1185 const G4double NKToNK2piCX = CrossSections::NKToNK2pi(particle1,particle2);
1186 const G4double totCX=CrossSections::total(particle1, particle2);
1187// assert(std::fabs(totCX-elasticCX-quasielasticCX-NKToNKpiCX-NKToNK2piCX)<0.1);
1188
1189 const G4double rChannel=Random::shoot() * totCX;
1190 if(elasticCX > rChannel){
1191// Elastic KN channel is chosen
1192 isElastic = true;
1193 INCL_DEBUG("KN interaction: elastic channel chosen" << '\n');
1194 return new NKElasticChannel(particle1, particle2);
1195 } else if(elasticCX + quasielasticCX > rChannel){
1196// Quasi-elastic KN channel is chosen
1197 isElastic = false; // true ??
1198 INCL_DEBUG("KN interaction: quasi-elastic channel chosen" << '\n');
1199 return new NKToNKChannel(particle1, particle2);
1200 } else if(elasticCX + quasielasticCX + NKToNKpiCX > rChannel){
1201// KN -> NKpi channel is chosen
1202 isElastic = false;
1203 INCL_DEBUG("KN interaction: NKpi channel chosen" << '\n');
1204 return new NKToNKpiChannel(particle1, particle2);
1205 } else if(elasticCX + quasielasticCX + NKToNKpiCX + NKToNK2piCX > rChannel){
1206// KN -> NK2pi channel is chosen
1207 isElastic = false;
1208 INCL_DEBUG("KN interaction: NK2pi channel chosen" << '\n');
1209 return new NKToNK2piChannel(particle1, particle2);
1210 } else {
1211 INCL_WARN("inconsistency within the KN Cross Sections (sum!=inelastic)" << '\n');
1212 if(NKToNK2piCX>0.) {
1213 INCL_WARN("Returning a NKToNK2pi channel" << '\n');
1214 isElastic = false;
1215 return new NKToNK2piChannel(particle1, particle2);
1216 } else if(NKToNKpiCX>0.) {
1217 INCL_WARN("Returning a NKToNKpi channel" << '\n');
1218 isElastic = false;
1219 return new NKToNKpiChannel(particle1, particle2);
1220 } else if(quasielasticCX>0.) {
1221 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1222 isElastic = false; // true ??
1223 return new NKToNKChannel(particle1, particle2);
1224 } else {
1225 INCL_WARN("Returning an elastic channel" << '\n');
1226 isElastic = true;
1227 return new NKElasticChannel(particle1, particle2);
1228 }
1229 }
1230 } else if ((particle1->isNucleon() && particle2->isAntiKaon()) || (particle2->isNucleon() && particle1->isAntiKaon())) {
1231//// KbN
1232 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1233 const G4double quasielasticCX = CrossSections::NKbToNKb(particle1,particle2);
1234 const G4double NKbToNKbpiCX = CrossSections::NKbToNKbpi(particle1,particle2);
1235 const G4double NKbToNKb2piCX = CrossSections::NKbToNKb2pi(particle1,particle2);
1236 const G4double NKbToLpiCX = CrossSections::NKbToLpi(particle1,particle2);
1237 const G4double NKbToL2piCX = CrossSections::NKbToL2pi(particle1,particle2);
1238 const G4double NKbToSpiCX = CrossSections::NKbToSpi(particle1,particle2);
1239 const G4double NKbToS2piCX = CrossSections::NKbToS2pi(particle1,particle2);
1240 const G4double totCX=CrossSections::total(particle1, particle2);
1241// assert(std::fabs(totCX-elasticCX-quasielasticCX-NKbToNKbpiCX-NKbToNKb2piCX-NKbToLpiCX-NKbToL2piCX-NKbToSpiCX-NKbToS2piCX)<0.1);
1242
1243 const G4double rChannel=Random::shoot() * totCX;
1244 if(elasticCX > rChannel){
1245// Elastic KbN channel is chosen
1246 isElastic = true;
1247 INCL_DEBUG("KbN interaction: elastic channel chosen" << '\n');
1248 return new NKbElasticChannel(particle1, particle2);
1249 } else if(elasticCX + quasielasticCX > rChannel){
1250// Quasi-elastic KbN channel is chosen
1251 isElastic = false; // true ??
1252 INCL_DEBUG("KbN interaction: quasi-elastic channel chosen" << '\n');
1253 return new NKbToNKbChannel(particle1, particle2);
1254 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX > rChannel){
1255// KbN -> NKbpi channel is chosen
1256 isElastic = false;
1257 INCL_DEBUG("KbN interaction: NKbpi channel chosen" << '\n');
1258 return new NKbToNKbpiChannel(particle1, particle2);
1259 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX > rChannel){
1260// KbN -> NKb2pi channel is chosen
1261 isElastic = false;
1262 INCL_DEBUG("KbN interaction: NKb2pi channel chosen" << '\n');
1263 return new NKbToNKb2piChannel(particle1, particle2);
1264 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX > rChannel){
1265// KbN -> Lpi channel is chosen
1266 isElastic = false;
1267 INCL_DEBUG("KbN interaction: Lpi channel chosen" << '\n');
1268 return new NKbToLpiChannel(particle1, particle2);
1269 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX > rChannel){
1270// KbN -> L2pi channel is chosen
1271 isElastic = false;
1272 INCL_DEBUG("KbN interaction: L2pi channel chosen" << '\n');
1273 return new NKbToL2piChannel(particle1, particle2);
1274 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX > rChannel){
1275// KbN -> Spi channel is chosen
1276 isElastic = false;
1277 INCL_DEBUG("KbN interaction: Spi channel chosen" << '\n');
1278 return new NKbToSpiChannel(particle1, particle2);
1279 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX + NKbToS2piCX > rChannel){
1280// KbN -> S2pi channel is chosen
1281 isElastic = false;
1282 INCL_DEBUG("KbN interaction: S2pi channel chosen" << '\n');
1283 return new NKbToS2piChannel(particle1, particle2);
1284 } else {
1285 INCL_WARN("inconsistency within the KbN Cross Sections (sum!=inelastic)" << '\n');
1286 if(NKbToS2piCX>0.) {
1287 INCL_WARN("Returning a NKbToS2pi channel" << '\n');
1288 isElastic = false;
1289 return new NKbToS2piChannel(particle1, particle2);
1290 } else if(NKbToSpiCX>0.) {
1291 INCL_WARN("Returning a NKbToSpi channel" << '\n');
1292 isElastic = false;
1293 return new NKbToSpiChannel(particle1, particle2);
1294 } else if(NKbToL2piCX>0.) {
1295 INCL_WARN("Returning a NKbToL2pi channel" << '\n');
1296 isElastic = false;
1297 return new NKbToL2piChannel(particle1, particle2);
1298 } else if(NKbToLpiCX>0.) {
1299 INCL_WARN("Returning a NKbToLpi channel" << '\n');
1300 isElastic = false;
1301 return new NKbToLpiChannel(particle1, particle2);
1302 } else if(NKbToNKb2piCX>0.) {
1303 INCL_WARN("Returning a NKbToNKb2pi channel" << '\n');
1304 isElastic = false;
1305 return new NKbToNKb2piChannel(particle1, particle2);
1306 } else if(NKbToNKbpiCX>0.) {
1307 INCL_WARN("Returning a NKbToNKbpi channel" << '\n');
1308 isElastic = false;
1309 return new NKbToNKbpiChannel(particle1, particle2);
1310 } else if(quasielasticCX>0.) {
1311 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1312 isElastic = false; // true ??
1313 return new NKbToNKbChannel(particle1, particle2);
1314 } else {
1315 INCL_WARN("Returning an elastic channel" << '\n');
1316 isElastic = true;
1317 return new NKbElasticChannel(particle1, particle2);
1318 }
1319 }
1320 } else if ((particle1->isNucleon() && particle2->isLambda()) || (particle2->isNucleon() && particle1->isLambda())) {
1321//// NLambda
1322 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1323 const G4double NLToNSCX = CrossSections::NLToNS(particle1,particle2);
1324 const G4double totCX=CrossSections::total(particle1, particle2);
1325// assert(std::fabs(totCX-elasticCX-NLToNSCX)<0.1);
1326
1327 const G4double rChannel=Random::shoot() * totCX;
1328 if(elasticCX > rChannel){
1329// Elastic NLambda channel is chosen
1330 isElastic = true;
1331 INCL_DEBUG("NLambda interaction: elastic channel chosen" << '\n');
1332 return new NYElasticChannel(particle1, particle2);
1333 } else if(elasticCX + NLToNSCX > rChannel){
1334// Quasi-elastic NLambda channel is chosen
1335 isElastic = false; // true ??
1336 INCL_DEBUG("NLambda interaction: quasi-elastic channel chosen" << '\n');
1337 return new NLToNSChannel(particle1, particle2);
1338 } else {
1339 INCL_WARN("inconsistency within the NLambda Cross Sections (sum!=inelastic)" << '\n');
1340 if(NLToNSCX>0.) {
1341 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1342 isElastic = false; // true ??
1343 return new NLToNSChannel(particle1, particle2);
1344 } else {
1345 INCL_WARN("Returning an elastic channel" << '\n');
1346 isElastic = true;
1347 return new NYElasticChannel(particle1, particle2);
1348 }
1349 }
1350 } else if ((particle1->isNucleon() && particle2->isSigma()) || (particle2->isNucleon() && particle1->isSigma())) {
1351//// NSigma
1352 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1353 const G4double NSToNLCX = CrossSections::NSToNL(particle1,particle2);
1354 const G4double NSToNSCX = CrossSections::NSToNS(particle1,particle2);
1355 const G4double totCX=CrossSections::total(particle1, particle2);
1356// assert(std::fabs(totCX-elasticCX-NSToNLCX-NSToNSCX)<0.1);
1357
1358 const G4double rChannel=Random::shoot() * totCX;
1359 if(elasticCX > rChannel){
1360// Elastic NSigma channel is chosen
1361 isElastic = true;
1362 INCL_DEBUG("NSigma interaction: elastic channel chosen" << '\n');
1363 return new NYElasticChannel(particle1, particle2);
1364 } else if(elasticCX + NSToNLCX > rChannel){
1365// NSigma -> NLambda channel is chosen
1366 isElastic = false; // true ??
1367 INCL_DEBUG("NSigma interaction: NLambda channel chosen" << '\n');
1368 return new NSToNLChannel(particle1, particle2);
1369 } else if(elasticCX + NSToNLCX + NSToNSCX > rChannel){
1370// NSigma -> NSigma quasi-elastic channel is chosen
1371 isElastic = false; // true ??
1372 INCL_DEBUG("NSigma interaction: NSigma quasi-elastic channel chosen" << '\n');
1373 return new NSToNSChannel(particle1, particle2);
1374 } else {
1375 INCL_WARN("inconsistency within the NSigma Cross Sections (sum!=inelastic)" << '\n');
1376 if(NSToNSCX>0.) {
1377 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1378 isElastic = false; // true ??
1379 return new NSToNSChannel(particle1, particle2);
1380 } else if(NSToNLCX>0.) {
1381 INCL_WARN("Returning a NLambda channel" << '\n');
1382 isElastic = false; // true ??
1383 return new NSToNLChannel(particle1, particle2);
1384 } else {
1385 INCL_WARN("Returning an elastic channel" << '\n');
1386 isElastic = true;
1387 return new NYElasticChannel(particle1, particle2);
1388 }
1389 }
1390 } else if ((particle1->isNucleon() && particle2->isAntiNucleon()) || (particle2->isNucleon() && particle1->isAntiNucleon())) {
1391//// NNbar
1392
1393 //Forcing annihilationfor emitInsideAntinucleon at the end of cascade && annihilation if E <= M + p_threhsold (from antideuteron in generateBinaryCollisionAvatar())
1394 /*const Particle *antinucleon;
1395 const Particle *nucleon;
1396
1397 if (particle1->isAntiNucleon()) {
1398 antinucleon = particle1;
1399 nucleon = particle2;
1400 }
1401 else {
1402 antinucleon = particle2;
1403 nucleon = particle1;
1404 }
1405 double Esquared = antinucleon->getEnergy() * antinucleon->getEnergy();
1406 double antinucleon_threshold=0;
1407 if(antinucleon->getType()==antiNeutron)
1408 antinucleon_threshold = theNucleus->getStore()->getConfig()->getnbAtrestThreshold();
1409 else if(antinucleon->getType() == antiProton)
1410 antinucleon_threshold = theNucleus->getStore()->getConfig()->getAtrestThreshold();
1411 else
1412 INCL_ERROR("neither antiproton nor antineutron");*/
1413 //double Sum_at_rest = antinucleon->getINCLMass() * antinucleon->getINCLMass() + antinucleon_threshold*antinucleon_threshold;
1414 // Force the annihilation when T < Threshold (at rest)
1415 //Config const *theConfig=theNucleus->getStore()->getConfig();
1416
1417 if ((theCrossSection == 9999.) || // XS=9999. means force annihilation
1418 ((particle1->getType()==antiProton && particle1->getKineticEnergy() <= theConfig->getAtrestThreshold()) ||
1419 (particle2->getType()==antiProton && particle2->getKineticEnergy() <= theConfig->getAtrestThreshold()) ||
1420 (particle1->getType()==antiNeutron && particle1->getKineticEnergy() <= particle1->getPotentialEnergy()) ||
1421 (particle2->getType()==antiNeutron && particle2->getKineticEnergy() <= particle2->getPotentialEnergy())) ||
1422 ((particle1->getType()==antiProton && particle1->getEnergy() <= particle1->getINCLMass()) ||
1423 (particle2->getType()==antiProton && particle2->getEnergy() <= particle2->getINCLMass()) ||
1424 (particle1->getType()==antiNeutron && particle1->getEnergy() <= particle1->getINCLMass()) ||
1425 (particle2->getType()==antiNeutron && particle2->getEnergy() <= particle2->getINCLMass())))
1426 {
1427 isElastic = false;
1428 AnnihilationType atype0;
1429 if((particle1->getType()==antiProton && particle2->getType()==Proton) || (particle2->getType()==antiProton && particle1->getType()==Proton)){
1430 atype0 = PTypeInFlight;
1431 }
1432 else if((particle1->getType()==antiProton && particle2->getType()==Neutron) || (particle2->getType()==antiProton && particle1->getType()==Neutron)){
1433 atype0 = NTypeInFlight;
1434 }
1435 else if((particle1->getType()==antiNeutron && particle2->getType()==Proton) || (particle2->getType()==antiNeutron && particle1->getType()==Proton)){
1436 atype0 = NbarPTypeInFlight;
1437 }
1438 else if((particle1->getType()==antiNeutron && particle2->getType()==Neutron) || (particle2->getType()==antiNeutron && particle1->getType()==Neutron)){
1439 atype0 = NbarNTypeInFlight;
1440 }
1441 else{
1442 atype0 = Def;
1443 INCL_ERROR("Annihilation type problem " << '\n');
1444 }
1445 theNucleus->setAType(atype0);
1446 return new NNbarToAnnihilationChannel(theNucleus, particle1, particle2);
1447 }
1448 //delete antinucleon;
1449 //delete nucleon;
1450
1451 // Usual interactions
1452 const G4double totCX = CrossSections::total(particle1, particle2);
1453 const G4double NNbElasticCX = CrossSections::NNbarElastic(particle1,particle2);
1454 const G4double NNbCEXCX = CrossSections::NNbarCEX(particle1,particle2);
1455 const G4double NNbToLLbCX = CrossSections::NNbarToLLbar(particle1,particle2);
1456 const G4double NNbToNNbpiCX = CrossSections::NNbarToNNbarpi(particle1,particle2);
1457 const G4double NNbToNNb2piCX = CrossSections::NNbarToNNbar2pi(particle1,particle2);
1458 const G4double NNbToNNb3piCX = CrossSections::NNbarToNNbar3pi(particle1,particle2);
1459 const G4double AnnihilationCX = CrossSections::NNbarToAnnihilation(particle1, particle2);
1460
1461// assert(std::fabs(totCX-NNbElasticCX-NNbCEXCX-NNbToLLbCX-NNbToNNbpiCX-NNbToNNb2piCX-NNbToNNb3piCX-AnnihilationCX)<0.1);
1462
1463 const G4double rChannel=Random::shoot() * totCX;
1464 if (NNbElasticCX > rChannel) {
1465 // NNbar (elastic) channel is chosen
1466 isElastic = true;
1467 //INCL_WARN("NNbar interaction: NNbarElastic channel chosen" << '\n');
1468 return new NNbarElasticChannel(particle1, particle2);
1469 } else if (NNbElasticCX + NNbCEXCX > rChannel) {
1470 // NNbar (CEX) channel is chosen
1471 isElastic = false; // may be charge-exchange also
1472 //INCL_WARN("NNbar interaction: NNbarCEX channel chosen" << '\n');
1473 return new NNbarCEXChannel(particle1, particle2);
1474 } else if (NNbElasticCX + NNbCEXCX + NNbToLLbCX > rChannel) {
1475 // NNbarToLLbar channel is chosen
1476 isElastic = false; // may be charge-exchange also
1477 //INCL_WARN("NNbar interaction: NNbarToLLbar channel chosen" << '\n');
1478 return new NNbarToLLbarChannel(particle1, particle2);
1479 } else if (NNbElasticCX + NNbCEXCX + NNbToLLbCX + NNbToNNbpiCX > rChannel) {
1480 // NNbar to NNbar pi channel is chosen
1481 isElastic = false;
1482 //INCL_WARN("NNbar interaction: NNbar pi channel chosen" << '\n');
1483 return new NNbarToNNbarpiChannel(particle1, particle2);
1484 } else if (NNbElasticCX + NNbCEXCX + NNbToLLbCX + NNbToNNbpiCX + NNbToNNb2piCX > rChannel) {
1485 // NNbar to NNbar 2pi channel is chosen
1486 isElastic = false;
1487 //INCL_WARN("NNbar interaction: NNbar 2pi channel chosen" << '\n');
1488 return new NNbarToNNbar2piChannel(particle1, particle2);
1489 } else if (NNbElasticCX + NNbCEXCX + NNbToLLbCX + NNbToNNbpiCX + NNbToNNb2piCX + NNbToNNb3piCX > rChannel) {
1490 // NNbar to NNbar 3pi channel is chosen
1491 isElastic = false;
1492 //INCL_WARN("NNbar interaction: NNbar 3pi channel chosen" << '\n');
1493 return new NNbarToNNbar3piChannel(particle1, particle2);
1494 } else if (NNbElasticCX + NNbCEXCX + NNbToLLbCX + NNbToNNbpiCX + NNbToNNb2piCX + NNbToNNb3piCX +AnnihilationCX > rChannel){
1495 // NNbar annihilation channel is chosen
1496 isElastic = false;
1497 AnnihilationType atype;
1498 if((particle1->getType()==antiProton && particle2->getType()==Proton) || (particle2->getType()==antiProton && particle1->getType()==Proton)){
1499 atype = PTypeInFlight;
1500 }
1501 else if((particle1->getType()==antiProton && particle2->getType()==Neutron) || (particle2->getType()==antiProton && particle1->getType()==Neutron)){
1502 atype = NTypeInFlight;
1503 }
1504 else if((particle1->getType()==antiNeutron && particle2->getType()==Proton) || (particle2->getType()==antiNeutron && particle1->getType()==Proton)){
1505 atype = NbarPTypeInFlight;
1506 }
1507 else if((particle1->getType()==antiNeutron && particle2->getType()==Neutron) || (particle2->getType()==antiNeutron && particle1->getType()==Neutron)){
1508 atype = NbarNTypeInFlight;
1509 }
1510 else{
1511 atype = Def;
1512 INCL_ERROR("Annihilation type problem " << '\n');
1513 }
1514 theNucleus->setAType(atype);
1515 return new NNbarToAnnihilationChannel(theNucleus, particle1, particle2);
1516 } else {
1517 INCL_WARN("Inconsistency within the NNbar Cross Sections (sum != inelastic)" << '\n');
1518 if (NNbToNNb3piCX > 0.0) {
1519 INCL_WARN("Returning an NNbar 3pi channel" << '\n');
1520 isElastic = false;
1521 return new NNbarToNNbar3piChannel(particle1, particle2);
1522 } else if (NNbToNNb2piCX > 0.0) {
1523 INCL_WARN("Returning an NNbar 2pi channel" << '\n');
1524 isElastic = false;
1525 return new NNbarToNNbar2piChannel(particle1, particle2);
1526 } else if (NNbToNNbpiCX > 0.0) {
1527 INCL_WARN("Returning an NNbar pi channel" << '\n');
1528 isElastic = false;
1529 return new NNbarToNNbarpiChannel(particle1, particle2);
1530 } else if (AnnihilationCX > 0.0) {
1531 INCL_WARN("Returning an NNbar annihilation channel" << '\n');
1532 isElastic = false;
1533 AnnihilationType atype;
1534 if((particle1->getType()==antiProton && particle2->getType()==Proton) || (particle2->getType()==antiProton && particle1->getType()==Proton)){
1535 atype = PTypeInFlight;
1536 }
1537 else if((particle1->getType()==antiProton && particle2->getType()==Neutron) || (particle2->getType()==antiProton && particle1->getType()==Neutron)){
1538 atype = NTypeInFlight;
1539 }
1540 else if((particle1->getType()==antiNeutron && particle2->getType()==Proton) || (particle2->getType()==antiNeutron && particle1->getType()==Proton)){
1541 atype = NbarPTypeInFlight;
1542 }
1543 else if((particle1->getType()==antiNeutron && particle2->getType()==Neutron) || (particle2->getType()==antiNeutron && particle1->getType()==Neutron)){
1544 atype = NbarNTypeInFlight;
1545 }
1546 else{
1547 atype = Def;
1548 INCL_ERROR("Annihilation type problem " << '\n');
1549 }
1550 theNucleus->setAType(atype);
1551 return new NNbarToAnnihilationChannel(theNucleus, particle1, particle2);
1552 } else if (NNbCEXCX > 0.0) {
1553 INCL_WARN("Returning an NNbar CEX channel" << '\n');
1554 isElastic = false;
1555 return new NNbarCEXChannel(particle1, particle2);
1556 } else if (NNbToLLbCX > 0.0) {
1557 INCL_WARN("Returning an NNbar LLbar channel" << '\n');
1558 isElastic = false;
1559 return new NNbarToLLbarChannel(particle1, particle2);
1560 } else {
1561 INCL_WARN("Elastic NNbar channel chosen" << '\n');
1562 isElastic = true;
1563 return new NNbarElasticChannel(particle1, particle2);
1564 }
1565 }
1566 }
1567
1568 else {
1569 INCL_DEBUG("BinaryCollisionAvatar can only handle nucleons (for the moment)."
1570 << '\n'
1571 << particle1->print()
1572 << '\n'
1573 << particle2->print()
1574 << '\n');
1575 InteractionAvatar::restoreParticles();
1576 return NULL;
1577 }
1578 }
INCL_ERROR
#define INCL_ERROR(x)
Definition G4INCLLogger.hh:225
INCL_WARN
#define INCL_WARN(x)
Definition G4INCLLogger.hh:232
INCL_DEBUG
#define INCL_DEBUG(x)
Definition G4INCLLogger.hh:240
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4INCL::InteractionAvatar::restoreParticles
void restoreParticles() const
Restore the state of both particles.
Definition G4INCLInteractionAvatar.cc:474
G4INCL::Particle::MergeVectorBias
static std::vector< G4int > MergeVectorBias(Particle const *const p1, Particle const *const p2)
Definition G4INCLParticle.cc:238
G4INCL::Particle::getBiasFromVector
static G4double getBiasFromVector(std::vector< G4int > VectorBias)
Definition G4INCLParticle.cc:226
G4INCL::CrossSections::NSToNL
G4double NSToNL(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:276
G4INCL::CrossSections::NNToNNKKb
G4double NNToNNKKb(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:203
G4INCL::CrossSections::elastic
G4double elastic(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:57
G4INCL::CrossSections::NpiToSK
G4double NpiToSK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:232
G4INCL::CrossSections::piNToOmegaN
G4double piNToOmegaN(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:89
G4INCL::CrossSections::NNbarCEX
G4double NNbarCEX(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:311
G4INCL::CrossSections::NKbToNKb2pi
G4double NKbToNKb2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:344
G4INCL::CrossSections::NNToNDeltaOmega
G4double NNToNDeltaOmega(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:162
G4INCL::CrossSections::NDeltaToDeltaLK
G4double NDeltaToDeltaLK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:217
G4INCL::CrossSections::NNToNSKpi
G4double NNToNSKpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:191
G4INCL::CrossSections::etaNToPiN
G4double etaNToPiN(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:97
G4INCL::CrossSections::NDeltaToNNKKb
G4double NDeltaToNNKKb(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:224
G4INCL::CrossSections::NKbToLpi
G4double NKbToLpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:304
G4INCL::CrossSections::NNToNLK2pi
G4double NNToNLK2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:195
G4INCL::CrossSections::etaNToPiPiN
G4double etaNToPiPiN(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:101
G4INCL::CrossSections::NKbToS2pi
G4double NKbToS2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:332
G4INCL::CrossSections::piNToEtaN
G4double piNToEtaN(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:85
G4INCL::CrossSections::NNbarToAnnihilation
G4double NNbarToAnnihilation(Particle const *const p1, Particle const *const p2)
Nucleon-AntiNucleon total annihilation cross sections.
Definition G4INCLCrossSections.cc:328
G4INCL::CrossSections::omegaNToPiN
G4double omegaNToPiN(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:113
G4INCL::CrossSections::NDeltaToNSK
G4double NDeltaToNSK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:214
G4INCL::CrossSections::NKbToSpi
G4double NKbToSpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:300
G4INCL::CrossSections::NNToNDelta
G4double NNToNDelta(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:69
G4INCL::CrossSections::NNToNSK
G4double NNToNSK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:183
G4INCL::CrossSections::NDeltaToNLK
G4double NDeltaToNLK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:211
G4INCL::CrossSections::NLToNS
G4double NLToNS(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:272
G4INCL::CrossSections::omegaNToLK
G4double omegaNToLK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:121
G4INCL::CrossSections::piNToDelta
G4double piNToDelta(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:77
G4INCL::CrossSections::NKbToNKb
G4double NKbToNKb(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:296
G4INCL::CrossSections::NNToNSK2pi
G4double NNToNSK2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:199
G4INCL::CrossSections::NNbarToNNbar3pi
G4double NNbarToNNbar3pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:324
G4INCL::CrossSections::NNToNLK
G4double NNToNLK(Particle const *const p1, Particle const *const p2)
Strange cross sections.
Definition G4INCLCrossSections.cc:179
G4INCL::CrossSections::NNToNNEtaExclu
G4double NNToNNEtaExclu(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:137
G4INCL::CrossSections::NNToNDeltaEta
G4double NNToNDeltaEta(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:145
G4INCL::CrossSections::NNToNNOmegaxPi
G4double NNToNNOmegaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:158
G4INCL::CrossSections::NNToNNEtaxPi
G4double NNToNNEtaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:141
G4INCL::CrossSections::NNbarToNNbarpi
G4double NNbarToNNbarpi(Particle const *const p1, Particle const *const p2)
Nucleon-AntiNucleon to Nucleon-AntiNucleon + pions cross sections.
Definition G4INCLCrossSections.cc:318
G4INCL::CrossSections::NSToNS
G4double NSToNS(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:280
G4INCL::CrossSections::NpiToNKKb
G4double NpiToNKKb(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:264
G4INCL::CrossSections::NpiToLK2pi
G4double NpiToLK2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:256
G4INCL::CrossSections::omegaNToPiPiN
G4double omegaNToPiPiN(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:117
G4INCL::CrossSections::NDeltaToDeltaSK
G4double NDeltaToDeltaSK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:220
G4INCL::CrossSections::NpiToMissingStrangeness
G4double NpiToMissingStrangeness(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:268
G4INCL::CrossSections::NNbarToNNbar2pi
G4double NNbarToNNbar2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:321
G4INCL::CrossSections::total
G4double total(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:61
G4INCL::CrossSections::NDeltaToNN
G4double NDeltaToNN(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:65
G4INCL::CrossSections::omegaNToSK
G4double omegaNToSK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:125
G4INCL::CrossSections::NpiToLKpi
G4double NpiToLKpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:248
G4INCL::CrossSections::NKbToNKbpi
G4double NKbToNKbpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:340
G4INCL::CrossSections::NNToNLKpi
G4double NNToNLKpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:187
G4INCL::CrossSections::NKbToL2pi
G4double NKbToL2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:336
G4INCL::CrossSections::NNToxPiNN
G4double NNToxPiNN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:73
G4INCL::CrossSections::NKToNK2pi
G4double NKToNK2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:292
G4INCL::CrossSections::NKToNKpi
G4double NKToNKpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:288
G4INCL::CrossSections::NNToNNOmegaExclu
G4double NNToNNOmegaExclu(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:154
G4INCL::CrossSections::NpiToSK2pi
G4double NpiToSK2pi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:260
G4INCL::CrossSections::NKToNK
G4double NKToNK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:284
G4INCL::CrossSections::NNToMissingStrangeness
G4double NNToMissingStrangeness(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:207
G4INCL::CrossSections::piNToxPiN
G4double piNToxPiN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:81
G4INCL::CrossSections::etaNToSK
G4double etaNToSK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:109
G4INCL::CrossSections::etaNToLK
G4double etaNToLK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:105
G4INCL::CrossSections::NNbarToLLbar
G4double NNbarToLLbar(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:314
G4INCL::CrossSections::NpiToLK
G4double NpiToLK(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:228
G4INCL::CrossSections::NNbarElastic
G4double NNbarElastic(Particle const *const p1, Particle const *const p2)
antiparticle cross sections
Definition G4INCLCrossSections.cc:308
G4INCL::CrossSections::NpiToSKpi
G4double NpiToSKpi(Particle const *const p1, Particle const *const p2)
Definition G4INCLCrossSections.cc:252
G4INCL::KinematicsUtils::squareTotalEnergyInCM
G4double squareTotalEnergyInCM(Particle const *const p1, Particle const *const p2)
Definition G4INCLKinematicsUtils.cc:125
G4INCL::Math::tenPi
const G4double tenPi
Definition G4INCLGlobals.hh:71
G4INCL::ParticleTable::getsrcPairConfig
G4bool getsrcPairConfig()
Get the configuration of src-pair correlations.
Definition G4INCLParticleTable.cc:1382
G4INCL::Random::shoot
G4double shoot()
Definition G4INCLRandom.cc:93
G4INCL::NbarPTypeInFlight
@ NbarPTypeInFlight
Definition G4INCLNucleus.hh:65
G4INCL::NbarNTypeInFlight
@ NbarNTypeInFlight
Definition G4INCLNucleus.hh:65
PoPI::ParticleList
std::vector< Base * > ParticleList
Definition PoPI.hpp:186

◆ getCutNN()

G4double G4INCL::BinaryCollisionAvatar::getCutNN ( )
inlinestatic

Definition at line 79 of file G4INCLBinaryCollisionAvatar.hh.

79{ return cutNN; }

◆ getCutNNSquared()

G4double G4INCL::BinaryCollisionAvatar::getCutNNSquared ( )
inlinestatic

Definition at line 81 of file G4INCLBinaryCollisionAvatar.hh.

81{ return cutNNSquared; }

Referenced by G4INCL::StandardPropagationModel::generateBinaryCollisionAvatar().

◆ getParticles()

ParticleList G4INCL::BinaryCollisionAvatar::getParticles ( ) const
inlinevirtual

Implements G4INCL::IAvatar.

Definition at line 62 of file G4INCLBinaryCollisionAvatar.hh.

62 {
63 ParticleList theParticleList;
64 theParticleList.push_back(particle1);
65 theParticleList.push_back(particle2);
66 return theParticleList;
67 };

◆ postInteraction()

void G4INCL::BinaryCollisionAvatar::postInteraction ( FinalState * fs)
virtual

Implements G4INCL::IAvatar.

Definition at line 1586 of file G4INCLBinaryCollisionAvatar.cc.

1586 {
1587 // Call the postInteraction method of the parent class
1588 // (provides Pauli blocking and enforces energy conservation)
1589 InteractionAvatar::postInteraction(fs);
1590
1591 switch(fs->getValidity()) {
1592 case PauliBlockedFS:
1593 theNucleus->getStore()->getBook().incrementBlockedCollisions();
1594 break;
1595 case NoEnergyConservationFS:
1596 case ParticleBelowFermiFS:
1597 case ParticleBelowZeroFS:
1598 break;
1599 case ValidFS:
1600 Book &theBook = theNucleus->getStore()->getBook();
1601 theBook.incrementAcceptedCollisions();
1602
1603 if(theBook.getAcceptedCollisions() == 1) {
1604 // Store time and cross section of the first collision
1605 G4double t = theBook.getCurrentTime();
1606 theBook.setFirstCollisionTime(t);
1607 theBook.setFirstCollisionXSec(oldXSec);
1608 // Increase the number of Kaon by 1
1609 if(isStrangeProduction) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()+1);
1610 // Store position and momentum of the spectator on the first
1611 // collision
1612 if((isParticle1Spectator && isParticle2Spectator) || (!isParticle1Spectator && !isParticle2Spectator)) {
1613 INCL_ERROR("First collision must be within a target spectator and a non-target spectator");
1614 }
1615 if(isParticle1Spectator) {
1616 theBook.setFirstCollisionSpectatorPosition(backupParticle1->getPosition().mag());
1617 theBook.setFirstCollisionSpectatorMomentum(backupParticle1->getMomentum().mag());
1618 } else {
1619 theBook.setFirstCollisionSpectatorPosition(backupParticle2->getPosition().mag());
1620 theBook.setFirstCollisionSpectatorMomentum(backupParticle2->getMomentum().mag());
1621 }
1622
1623 // Store the elasticity of the first collision
1624 theBook.setFirstCollisionIsElastic(isElastic);
1625 }
1626 }
1627 return;
1628 }
G4INCL::InteractionAvatar::backupParticle2
static G4ThreadLocal Particle * backupParticle2
Definition G4INCLInteractionAvatar.hh:123
G4INCL::InteractionAvatar::backupParticle1
static G4ThreadLocal Particle * backupParticle1
Definition G4INCLInteractionAvatar.hh:123
G4INCL::NoEnergyConservationFS
@ NoEnergyConservationFS
Definition G4INCLFinalState.hh:56

◆ preInteraction()

void G4INCL::BinaryCollisionAvatar::preInteraction ( )
virtual

Implements G4INCL::IAvatar.

Definition at line 1580 of file G4INCLBinaryCollisionAvatar.cc.

1580 {
1581 isParticle1Spectator = particle1->isTargetSpectator();
1582 isParticle2Spectator = particle2->isTargetSpectator();
1583 InteractionAvatar::preInteraction();
1584 }

◆ setBias()

void G4INCL::BinaryCollisionAvatar::setBias ( const G4double b)
inlinestatic

Set the global bias factor.

Definition at line 87 of file G4INCLBinaryCollisionAvatar.hh.

87{ bias=b; }

Referenced by G4INCL::INCL::INCL().

◆ setCutNN()

void G4INCL::BinaryCollisionAvatar::setCutNN ( const G4double c)
inlinestatic

Definition at line 74 of file G4INCLBinaryCollisionAvatar.hh.

74 {
75 cutNN = c;
76 cutNNSquared = cutNN*cutNN;
77 }

Referenced by G4INCL::INCL::INCL().

The documentation for this class was generated from the following files: