Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | List of all members
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 ()
 
- Public Member Functions inherited from G4INCL::IAvatar
 IAvatar ()
 
 IAvatar (G4double time)
 
virtual ~IAvatar ()
 
virtual G4INCL::IChannel * getChannel ()=0
 
FinalState * getFinalState ()
 
void fillFinalState (FinalState *fs)
 
virtual void preInteraction ()=0
 
virtual void postInteraction (FinalState *)=0
 
G4double getTime () const
 
virtual ParticleList getParticles () const =0
 
virtual std::string dump () const =0
 
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.
 

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
virtual G4INCL::IChannel * getChannel ()=0
 
G4bool bringParticleInside (Particle *const p)
 
void preInteractionLocalEnergy (Particle *const p)
 Apply local-energy transformation, if appropriate.
 
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.
 
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
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 56 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 124 of file G4INCLBinaryCollisionAvatar.cc.

126 : InteractionAvatar(time, n, p1, p2), theCrossSection(crossSection),
127 isParticle1Spectator(false),
128 isParticle2Spectator(false),
129 isElastic(false),
130 isStrangeProduction(false)
131 {
132 setType(CollisionAvatarType);
133 }
G4INCL::IAvatar::setType
void setType(AvatarType t)
Definition: G4INCLIAvatar.hh:93
G4INCL::InteractionAvatar::InteractionAvatar
InteractionAvatar(G4double, G4INCL::Nucleus *, G4INCL::Particle *)
Definition: G4INCLInteractionAvatar.cc:65
G4INCL::CollisionAvatarType
@ CollisionAvatarType
Definition: G4INCLIAvatar.hh:67

◆ ~BinaryCollisionAvatar()

G4INCL::BinaryCollisionAvatar::~BinaryCollisionAvatar ( )
virtual

Definition at line 135 of file G4INCLBinaryCollisionAvatar.cc.

135 {
136 }

Member Function Documentation

◆ dump()

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

Implements G4INCL::IAvatar.

Definition at line 1300 of file G4INCLBinaryCollisionAvatar.cc.

1300 {
1301 std::stringstream ss;
1302 ss << "(avatar " << theTime <<" 'nn-collision" << '\n'
1303 << "(list " << '\n'
1304 << particle1->dump()
1305 << particle2->dump()
1306 << "))" << '\n';
1307 return ss.str();
1308 }
G4INCL::Particle::dump
std::string dump() const
Definition: G4INCLParticle.hh:972

◆ getBias()

static G4double G4INCL::BinaryCollisionAvatar::getBias ( )
inlinestatic

Get the global bias factor.

Definition at line 83 of file G4INCLBinaryCollisionAvatar.hh.

83{ 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 138 of file G4INCLBinaryCollisionAvatar.cc.

138 {
139 // We already check cutNN at avatar creation time, but we have to check it
140 // again here. For composite projectiles, we might have created independent
141 // avatars with no cutNN before any collision took place.
142 if(particle1->isNucleon()
143 && particle2->isNucleon()
144 && theNucleus->getStore()->getBook().getAcceptedCollisions()!=0) {
145 const G4double energyCM2 = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
146 // Below a certain cut value we don't do anything:
147 if(energyCM2 < cutNNSquared) {
148 INCL_DEBUG("CM energy = sqrt(" << energyCM2 << ") MeV < std::sqrt(" << cutNNSquared
149 << ") MeV = cutNN" << "; returning a NULL channel" << '\n');
150 InteractionAvatar::restoreParticles();
151 return NULL;
152 }
153 }
154
155 /** Check again the distance of approach. In order for the avatar to be
156 * realised, we have to perform a check in the CM system. We define a
157 * distance four-vector as
158 * \f[ (0, \Delta\vec{x}), \f]
159 * where \f$\Delta\vec{x}\f$ is the distance vector of the particles at
160 * their minimum distance of approach (i.e. at the avatar time). By
161 * boosting this four-vector to the CM frame of the two particles and we
162 * obtain a new four vector
163 * \f[ (\Delta t', \Delta\vec{x}'), \f]
164 * with a non-zero time component (the collision happens simultaneously for
165 * the two particles in the lab system, but not in the CM system). In order
166 * for the avatar to be realised, we require that
167 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/\pi}.\f]
168 * Note that \f$|\Delta\vec{x}'|\leq|\Delta\vec{x}|\f$; thus, the condition
169 * above is more restrictive than the check that we perform in
170 * G4INCL::Propagation::StandardPropagationModel::generateBinaryCollisionAvatar.
171 * In other words, the avatar generation cannot miss any physical collision
172 * avatars.
173 */
174 ThreeVector minimumDistance = particle1->getPosition();
175 minimumDistance -= particle2->getPosition();
176 const G4double betaDotX = boostVector.dot(minimumDistance);
177 const G4double minDist = Math::tenPi*(minimumDistance.mag2() + betaDotX*betaDotX / (1.-boostVector.mag2()));
178 if(minDist > theCrossSection) {
179 INCL_DEBUG("CM distance of approach is too small: " << minDist << ">" <<
180 theCrossSection <<"; returning a NULL channel" << '\n');
181 InteractionAvatar::restoreParticles();
182 return NULL;
183 }
184
185 /** Bias apply for this reaction in order to get the same
186 * ParticleBias for all stange particles.
187 * Can be reduced after because of the safeguard.
188 */
189 G4double bias_apply = 1.;
190 if(bias != 1.) bias_apply = Particle::getBiasFromVector(Particle::MergeVectorBias(particle1,particle2)) * bias;
191
192//// NN
193 if(particle1->isNucleon() && particle2->isNucleon()) {
194
195 G4double NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
196 G4double NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
197 G4double NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
198 G4double NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
199 G4double NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
200 G4double NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
201 G4double NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
202 G4double NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
203
204 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NNToNDelta(particle1, particle2) + CrossSections::NNToxPiNN(1,particle1, particle2)
205 + CrossSections::NNToxPiNN(2,particle1, particle2) + CrossSections::NNToxPiNN(3,particle1, particle2) + CrossSections::NNToxPiNN(4,particle1, particle2)
206 + CrossSections::NNToNNEtaExclu(particle1, particle2) + CrossSections::NNToNDeltaEta(particle1, particle2) + CrossSections::NNToNNEtaxPi(1,particle1, particle2)
207 + CrossSections::NNToNNEtaxPi(2,particle1, particle2) + CrossSections::NNToNNEtaxPi(3,particle1, particle2) + CrossSections::NNToNNEtaxPi(4,particle1, particle2)
208 + CrossSections::NNToNNOmegaExclu(particle1, particle2) + CrossSections::NNToNDeltaOmega(particle1, particle2) + CrossSections::NNToNNOmegaxPi(1,particle1, particle2)
209 + CrossSections::NNToNNOmegaxPi(2,particle1, particle2) + CrossSections::NNToNNOmegaxPi(3,particle1, particle2) + CrossSections::NNToNNOmegaxPi(4,particle1, particle2);
210 const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + NLKpiProductionCX + NSKpiProductionCX + NLK2piProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX)/bias_apply;
211
212 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
213
214 if(counterweight < 0.5) {
215 counterweight = 0.5;
216 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
217 NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
218 NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
219 NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
220 NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
221 NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
222 NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
223 NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
224 NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
225 }
226
227
228 const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
229 const G4double deltaProductionCX = CrossSections::NNToNDelta(particle1, particle2)*counterweight;
230 const G4double onePiProductionCX = CrossSections::NNToxPiNN(1,particle1, particle2)*counterweight;
231 const G4double twoPiProductionCX = CrossSections::NNToxPiNN(2,particle1, particle2)*counterweight;
232 const G4double threePiProductionCX = CrossSections::NNToxPiNN(3,particle1, particle2)*counterweight;
233 const G4double fourPiProductionCX = CrossSections::NNToxPiNN(4,particle1, particle2)*counterweight;
234
235 const G4double etaProductionCX = CrossSections::NNToNNEtaExclu(particle1, particle2)*counterweight;
236 const G4double etadeltaProductionCX = CrossSections::NNToNDeltaEta(particle1, particle2)*counterweight;
237 const G4double etaonePiProductionCX = CrossSections::NNToNNEtaxPi(1,particle1, particle2)*counterweight;
238 const G4double etatwoPiProductionCX = CrossSections::NNToNNEtaxPi(2,particle1, particle2)*counterweight;
239 const G4double etathreePiProductionCX = CrossSections::NNToNNEtaxPi(3,particle1, particle2)*counterweight;
240 const G4double etafourPiProductionCX = CrossSections::NNToNNEtaxPi(4,particle1, particle2)*counterweight;
241 const G4double omegaProductionCX = CrossSections::NNToNNOmegaExclu(particle1, particle2)*counterweight;
242 const G4double omegadeltaProductionCX = CrossSections::NNToNDeltaOmega(particle1, particle2)*counterweight;
243 const G4double omegaonePiProductionCX = CrossSections::NNToNNOmegaxPi(1,particle1, particle2)*counterweight;
244 const G4double omegatwoPiProductionCX = CrossSections::NNToNNOmegaxPi(2,particle1, particle2)*counterweight;
245 const G4double omegathreePiProductionCX = CrossSections::NNToNNOmegaxPi(3,particle1, particle2)*counterweight;
246 const G4double omegafourPiProductionCX = CrossSections::NNToNNOmegaxPi(4,particle1, particle2)*counterweight;
247
248 const G4double totCX=CrossSections::total(particle1, particle2);
249
250// 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);
251
252 const G4double rChannel=Random::shoot() * totCX;
253
254 if(elasticCX > rChannel) {
255// Elastic NN channel
256 isElastic = true;
257 INCL_DEBUG("NN interaction: elastic channel chosen" << '\n');
258 weight = counterweight;
259 return new ElasticChannel(particle1, particle2);
260 } else if((elasticCX + deltaProductionCX) > rChannel) {
261 isElastic = false;
262// NN -> N Delta channel is chosen
263 INCL_DEBUG("NN interaction: Delta channel chosen" << '\n');
264 weight = counterweight;
265 return new DeltaProductionChannel(particle1, particle2);
266 } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
267 isElastic = false;
268// NN -> PiNN channel is chosen
269 INCL_DEBUG("NN interaction: one Pion channel chosen" << '\n');
270 weight = counterweight;
271 return new NNToMultiPionsChannel(1,particle1, particle2);
272 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
273 isElastic = false;
274// NN -> 2PiNN channel is chosen
275 INCL_DEBUG("NN interaction: two Pions channel chosen" << '\n');
276 weight = counterweight;
277 return new NNToMultiPionsChannel(2,particle1, particle2);
278 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
279 isElastic = false;
280// NN -> 3PiNN channel is chosen
281 INCL_DEBUG("NN interaction: three Pions channel chosen" << '\n');
282 weight = counterweight;
283 return new NNToMultiPionsChannel(3,particle1, particle2);
284 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX > rChannel) {
285 isElastic = false;
286// NN -> 4PiNN channel is chosen
287 INCL_DEBUG("NN interaction: four Pions channel chosen" << '\n');
288 weight = counterweight;
289 return new NNToMultiPionsChannel(4,particle1, particle2);
290 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
291 + etaProductionCX > rChannel) {
292 isElastic = false;
293// NN -> NNEta channel is chosen
294 INCL_DEBUG("NN interaction: Eta channel chosen" << '\n');
295 weight = counterweight;
296 return new NNToNNEtaChannel(particle1, particle2);
297 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
298 + etaProductionCX + etadeltaProductionCX > rChannel) {
299 isElastic = false;
300// NN -> N Delta Eta channel is chosen
301 INCL_DEBUG("NN interaction: Delta Eta channel chosen" << '\n');
302 weight = counterweight;
303 return new NDeltaEtaProductionChannel(particle1, particle2);
304 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
305 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX > rChannel) {
306 isElastic = false;
307// NN -> EtaPiNN channel is chosen
308 INCL_DEBUG("NN interaction: Eta + one Pion channel chosen" << '\n');
309 weight = counterweight;
310 return new NNEtaToMultiPionsChannel(1,particle1, particle2);
311 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
312 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX > rChannel) {
313 isElastic = false;
314// NN -> Eta2PiNN channel is chosen
315 INCL_DEBUG("NN interaction: Eta + two Pions channel chosen" << '\n');
316 weight = counterweight;
317 return new NNEtaToMultiPionsChannel(2,particle1, particle2);
318 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
319 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX > rChannel) {
320 isElastic = false;
321// NN -> Eta3PiNN channel is chosen
322 INCL_DEBUG("NN interaction: Eta + three Pions channel chosen" << '\n');
323 weight = counterweight;
324 return new NNEtaToMultiPionsChannel(3,particle1, particle2);
325 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
326 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX > rChannel) {
327 isElastic = false;
328// NN -> Eta4PiNN channel is chosen
329 INCL_DEBUG("NN interaction: Eta + four Pions channel chosen" << '\n');
330 weight = counterweight;
331 return new NNEtaToMultiPionsChannel(4,particle1, particle2);
332 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
333 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
334 + omegaProductionCX > rChannel) {
335 isElastic = false;
336// NN -> NNOmega channel is chosen
337 INCL_DEBUG("NN interaction: Omega channel chosen" << '\n');
338 weight = counterweight;
339 return new NNToNNOmegaChannel(particle1, particle2);
340 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
341 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
342 + omegaProductionCX + omegadeltaProductionCX > rChannel) {
343 isElastic = false;
344// NN -> N Delta Omega channel is chosen
345 INCL_DEBUG("NN interaction: Delta Omega channel chosen" << '\n');
346 weight = counterweight;
347 return new NDeltaOmegaProductionChannel(particle1, particle2);
348 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
349 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
350 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX > rChannel) {
351 isElastic = false;
352// NN -> OmegaPiNN channel is chosen
353 INCL_DEBUG("NN interaction: Omega + one Pion channel chosen" << '\n');
354 weight = counterweight;
355 return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
356 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
357 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
358 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX > rChannel) {
359 isElastic = false;
360// NN -> Omega2PiNN channel is chosen
361 INCL_DEBUG("NN interaction: Omega + two Pions channel chosen" << '\n');
362 weight = counterweight;
363 return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
364 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
365 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
366 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX > rChannel) {
367 isElastic = false;
368// NN -> Omega3PiNN channel is chosen
369 INCL_DEBUG("NN interaction: Omega + three Pions channel chosen" << '\n');
370 weight = counterweight;
371 return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
372 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
373 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
374 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX > rChannel) {
375 isElastic = false;
376// NN -> Omega4PiNN channel is chosen
377 INCL_DEBUG("NN interaction: Omega + four Pions channel chosen" << '\n');
378 weight = counterweight;
379 return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
380 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
381 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
382 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
383 + NLKProductionCX > rChannel) {
384 isElastic = false;
385 isStrangeProduction = true;
386// NN -> NLK channel is chosen
387 INCL_DEBUG("NN interaction: NLK channel chosen" << '\n');
388 weight = bias_apply;
389 return new NNToNLKChannel(particle1, particle2);
390 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
391 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
392 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
393 + NLKProductionCX + NLKpiProductionCX > rChannel) {
394 isElastic = false;
395 isStrangeProduction = true;
396// NN -> NLKpi channel is chosen
397 INCL_DEBUG("NN interaction: NLKpi channel chosen" << '\n');
398 weight = bias_apply;
399 return new NNToNLKpiChannel(particle1, particle2);
400 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
401 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
402 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
403 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX > rChannel) {
404 isElastic = false;
405 isStrangeProduction = true;
406// NN -> NLK2pi channel is chosen
407 INCL_DEBUG("NN interaction: NLK2pi channel chosen" << '\n');
408 weight = bias_apply;
409 return new NNToNLK2piChannel(particle1, particle2);
410 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
411 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
412 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
413 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX > rChannel) {
414 isElastic = false;
415 isStrangeProduction = true;
416// NN -> NSK channel is chosen
417 INCL_DEBUG("NN interaction: NSK channel chosen" << '\n');
418 weight = bias_apply;
419 return new NNToNSKChannel(particle1, particle2);
420 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
421 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
422 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
423 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX > rChannel) {
424 isElastic = false;
425 isStrangeProduction = true;
426// NN -> NSKpi channel is chosen
427 INCL_DEBUG("NN interaction: NSKpi channel chosen" << '\n');
428 weight = bias_apply;
429 return new NNToNSKpiChannel(particle1, particle2);
430 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
431 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
432 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
433 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX > rChannel) {
434 isElastic = false;
435 isStrangeProduction = true;
436// NN -> NSK2pi channel is chosen
437 INCL_DEBUG("NN interaction: NSK2pi channel chosen" << '\n');
438 weight = bias_apply;
439 return new NNToNSK2piChannel(particle1, particle2);
440 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
441 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
442 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
443 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX > rChannel) {
444 isElastic = false;
445 isStrangeProduction = true;
446// NN -> NNKKb channel is chosen
447 INCL_DEBUG("NN interaction: NNKKb channel chosen" << '\n');
448 weight = bias_apply;
449 return new NNToNNKKbChannel(particle1, particle2);
450 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
451 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
452 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
453 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX> rChannel) {
454 isElastic = false;
455 isStrangeProduction = true;
456// NN -> Missing Strangeness channel is chosen
457 INCL_DEBUG("NN interaction: Missing Strangeness channel chosen" << '\n');
458 weight = bias_apply;
459 return new NNToMissingStrangenessChannel(particle1, particle2);
460 } else {
461 INCL_WARN("inconsistency within the NN Cross Sections (sum!=inelastic)" << '\n');
462 if(NNMissingCX>0.) {
463 INCL_WARN("Returning an Missing Strangeness channel" << '\n');
464 weight = bias_apply;
465 isElastic = false;
466 isStrangeProduction = true;
467 return new NNToNNKKbChannel(particle1, particle2);
468 } else if(NNKKbProductionCX>0.) {
469 INCL_WARN("Returning an NNKKb channel" << '\n');
470 weight = bias_apply;
471 isElastic = false;
472 isStrangeProduction = true;
473 return new NNToNNKKbChannel(particle1, particle2);
474 } else if(NSK2piProductionCX>0.) {
475 INCL_WARN("Returning an NSK2pi channel" << '\n');
476 weight = bias_apply;
477 isElastic = false;
478 isStrangeProduction = true;
479 return new NNToNSK2piChannel(particle1, particle2);
480 } else if(NSKpiProductionCX>0.) {
481 INCL_WARN("Returning an NSKpi channel" << '\n');
482 weight = bias_apply;
483 isElastic = false;
484 isStrangeProduction = true;
485 return new NNToNSKpiChannel(particle1, particle2);
486 } else if(NSKProductionCX>0.) {
487 INCL_WARN("Returning an NSK channel" << '\n');
488 weight = bias_apply;
489 isElastic = false;
490 isStrangeProduction = true;
491 return new NNToNSKChannel(particle1, particle2);
492 } else if(NLK2piProductionCX>0.) {
493 INCL_WARN("Returning an NLK2pi channel" << '\n');
494 weight = bias_apply;
495 isElastic = false;
496 isStrangeProduction = true;
497 return new NNToNLK2piChannel(particle1, particle2);
498 } else if(NLKpiProductionCX>0.) {
499 INCL_WARN("Returning an NLKpi channel" << '\n');
500 weight = bias_apply;
501 isElastic = false;
502 isStrangeProduction = true;
503 return new NNToNLKpiChannel(particle1, particle2);
504 } else if(NLKProductionCX>0.) {
505 INCL_WARN("Returning an NLK channel" << '\n');
506 weight = bias_apply;
507 isElastic = false;
508 isStrangeProduction = true;
509 return new NNToNLKChannel(particle1, particle2);
510 } else if(omegafourPiProductionCX>0.) {
511 INCL_WARN("Returning an Omega + four Pions channel" << '\n');
512 weight = counterweight;
513 isElastic = false;
514 return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
515 } else if(omegathreePiProductionCX>0.) {
516 INCL_WARN("Returning an Omega + three Pions channel" << '\n');
517 weight = counterweight;
518 isElastic = false;
519 return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
520 } else if(omegatwoPiProductionCX>0.) {
521 INCL_WARN("Returning an Omega + two Pions channel" << '\n');
522 weight = counterweight;
523 isElastic = false;
524 return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
525 } else if(omegaonePiProductionCX>0.) {
526 INCL_WARN("Returning an Omega + one Pion channel" << '\n');
527 weight = counterweight;
528 isElastic = false;
529 return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
530 } else if(omegadeltaProductionCX>0.) {
531 INCL_WARN("Returning an Omega + Delta channel" << '\n');
532 weight = counterweight;
533 isElastic = false;
534 return new NDeltaOmegaProductionChannel(particle1, particle2);
535 } else if(omegaProductionCX>0.) {
536 INCL_WARN("Returning an Omega channel" << '\n');
537 weight = counterweight;
538 isElastic = false;
539 return new NNToNNOmegaChannel(particle1, particle2);
540 } else if(etafourPiProductionCX>0.) {
541 INCL_WARN("Returning an Eta + four Pions channel" << '\n');
542 weight = counterweight;
543 isElastic = false;
544 return new NNEtaToMultiPionsChannel(4,particle1, particle2);
545 } else if(etathreePiProductionCX>0.) {
546 INCL_WARN("Returning an Eta + threev channel" << '\n');
547 weight = counterweight;
548 isElastic = false;
549 return new NNEtaToMultiPionsChannel(3,particle1, particle2);
550 } else if(etatwoPiProductionCX>0.) {
551 INCL_WARN("Returning an Eta + two Pions channel" << '\n');
552 weight = counterweight;
553 isElastic = false;
554 return new NNEtaToMultiPionsChannel(2,particle1, particle2);
555 } else if(etaonePiProductionCX>0.) {
556 INCL_WARN("Returning an Eta + one Pion channel" << '\n');
557 weight = counterweight;
558 isElastic = false;
559 return new NNEtaToMultiPionsChannel(1,particle1, particle2);
560 } else if(etadeltaProductionCX>0.) {
561 INCL_WARN("Returning an Eta + Delta channel" << '\n');
562 weight = counterweight;
563 isElastic = false;
564 return new NDeltaEtaProductionChannel(particle1, particle2);
565 } else if(etaProductionCX>0.) {
566 INCL_WARN("Returning an Eta channel" << '\n');
567 weight = counterweight;
568 isElastic = false;
569 return new NNToNNEtaChannel(particle1, particle2);
570 } else if(fourPiProductionCX>0.) {
571 INCL_WARN("Returning a 4pi channel" << '\n');
572 weight = counterweight;
573 isElastic = false;
574 return new NNToMultiPionsChannel(4,particle1, particle2);
575 } else if(threePiProductionCX>0.) {
576 INCL_WARN("Returning a 3pi channel" << '\n');
577 weight = counterweight;
578 isElastic = false;
579 return new NNToMultiPionsChannel(3,particle1, particle2);
580 } else if(twoPiProductionCX>0.) {
581 INCL_WARN("Returning a 2pi channel" << '\n');
582 weight = counterweight;
583 isElastic = false;
584 return new NNToMultiPionsChannel(2,particle1, particle2);
585 } else if(onePiProductionCX>0.) {
586 INCL_WARN("Returning a 1pi channel" << '\n');
587 weight = counterweight;
588 isElastic = false;
589 return new NNToMultiPionsChannel(1,particle1, particle2);
590 } else if(deltaProductionCX>0.) {
591 INCL_WARN("Returning a delta-production channel" << '\n');
592 weight = counterweight;
593 isElastic = false;
594 return new DeltaProductionChannel(particle1, particle2);
595 } else {
596 INCL_WARN("Returning an elastic channel" << '\n');
597 weight = counterweight;
598 isElastic = true;
599 return new ElasticChannel(particle1, particle2);
600 }
601 }
602
603//// NDelta
604 }
605 else if((particle1->isNucleon() && particle2->isDelta()) ||
606 (particle1->isDelta() && particle2->isNucleon())) {
607
608 G4double NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
609 G4double NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
610 G4double DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
611 G4double DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
612 G4double NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
613
614 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NDeltaToNN(particle1, particle2);
615 const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX)/bias_apply;
616
617 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
618
619 if(counterweight < 0.5){
620 counterweight = 0.5;
621 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
622
623 NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
624 NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
625 DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
626 DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
627 NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
628 }
629
630 G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
631 G4double recombinationCX = CrossSections::NDeltaToNN(particle1, particle2)*counterweight;
632
633 const G4double rChannel=Random::shoot() * (StrangenessProdCX + UnStrangeProdCX);
634
635 if(elasticCX > rChannel) {
636 isElastic = true;
637// Elastic N Delta channel
638 INCL_DEBUG("NDelta interaction: elastic channel chosen" << '\n');
639 weight = counterweight;
640 return new ElasticChannel(particle1, particle2);
641 } else if (elasticCX + recombinationCX > rChannel){
642 isElastic = false;
643// Recombination
644// NDelta -> NN channel is chosen
645 INCL_DEBUG("NDelta interaction: recombination channel chosen" << '\n');
646 weight = counterweight;
647 return new RecombinationChannel(particle1, particle2);
648 } else if (elasticCX + recombinationCX + NLKProductionCX > rChannel){
649 isElastic = false;
650 isStrangeProduction = true;
651// NDelta -> NLK channel is chosen
652 INCL_DEBUG("NDelta interaction: NLK channel chosen" << '\n');
653 weight = bias_apply;
654 return new NDeltaToNLKChannel(particle1, particle2);
655 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX > rChannel){
656 isElastic = false;
657 isStrangeProduction = true;
658// NDelta -> NSK channel is chosen
659 INCL_DEBUG("NDelta interaction: NSK channel chosen" << '\n');
660 weight = bias_apply;
661 return new NDeltaToNSKChannel(particle1, particle2);
662 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX > rChannel){
663 isElastic = false;
664 isStrangeProduction = true;
665// NDelta -> DeltaLK channel is chosen
666 INCL_DEBUG("NDelta interaction: DeltaLK channel chosen" << '\n');
667 weight = bias_apply;
668 return new NDeltaToDeltaLKChannel(particle1, particle2);
669 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX > rChannel){
670 isElastic = false;
671 isStrangeProduction = true;
672// NDelta -> DeltaSK channel is chosen
673 INCL_DEBUG("NDelta interaction: DeltaSK channel chosen" << '\n');
674 weight = bias_apply;
675 return new NDeltaToDeltaSKChannel(particle1, particle2);
676 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX > rChannel){
677 isElastic = false;
678 isStrangeProduction = true;
679// NDelta -> NNKKb channel is chosen
680 INCL_DEBUG("NDelta interaction: NNKKb channel chosen" << '\n');
681 weight = bias_apply;
682 return new NDeltaToNNKKbChannel(particle1, particle2);
683 }
684 else{
685 INCL_ERROR("rChannel > (StrangenessProdCX + UnStrangeProdCX) in NDelta interaction: return an elastic channel" << '\n');
686 weight = counterweight;
687 isElastic = true;
688 return new ElasticChannel(particle1, particle2);
689 }
690
691//// DeltaDelta
692 } else if(particle1->isDelta() && particle2->isDelta()) {
693 isElastic = true;
694 INCL_DEBUG("DeltaDelta interaction: elastic channel chosen" << '\n');
695 return new ElasticChannel(particle1, particle2);
696
697//// PiN
698 } else if(isPiN) {
699
700 G4double LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
701 G4double SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
702 G4double LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
703 G4double SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
704 G4double LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
705 G4double SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
706 G4double NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
707 G4double MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
708
709 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::piNToDelta(particle1, particle2)
710 + CrossSections::piNToxPiN(2,particle1, particle2) + CrossSections::piNToxPiN(3,particle1, particle2) + CrossSections::piNToxPiN(4,particle1, particle2)
711 + CrossSections::piNToEtaN(particle1, particle2) + CrossSections::piNToOmegaN(particle1, particle2);
712 const G4double StrangenessProdCX = (LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX)/bias_apply;
713
714 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
715
716 if(counterweight < 0.5) {
717 counterweight = 0.5;
718 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
719 LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
720 SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
721 LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
722 SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
723 LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
724 SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
725 NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
726 MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
727 }
728
729
730 const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
731 const G4double deltaProductionCX = CrossSections::piNToDelta(particle1, particle2)*counterweight;
732 const G4double onePiProductionCX = CrossSections::piNToxPiN(2,particle1, particle2)*counterweight;
733 const G4double twoPiProductionCX = CrossSections::piNToxPiN(3,particle1, particle2)*counterweight;
734 const G4double threePiProductionCX = CrossSections::piNToxPiN(4,particle1, particle2)*counterweight;
735 const G4double etaProductionCX = CrossSections::piNToEtaN(particle1, particle2)*counterweight;
736 const G4double omegaProductionCX = CrossSections::piNToOmegaN(particle1, particle2)*counterweight;
737
738 const G4double totCX=CrossSections::total(particle1, particle2);
739
740// assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-etaProductionCX-omegaProductionCX-LKProdCX-SKProdCX-LKpiProdCX-SKpiProdCX-LK2piProdCX-SK2piProdCX-NKKbProdCX-MissingCX) < 0.15);
741
742 const G4double rChannel=Random::shoot() * totCX;
743
744 if(elasticCX > rChannel) {
745 isElastic = true;
746// Elastic PiN channel
747 INCL_DEBUG("PiN interaction: elastic channel chosen" << '\n');
748 weight = counterweight;
749 return new PiNElasticChannel(particle1, particle2);
750 } else if(elasticCX + deltaProductionCX > rChannel) {
751 isElastic = false;
752// PiN -> Delta channel is chosen
753 INCL_DEBUG("PiN interaction: Delta channel chosen" << '\n');
754 weight = counterweight;
755 return new PiNToDeltaChannel(particle1, particle2);
756 } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
757 isElastic = false;
758// PiN -> PiNPi channel is chosen
759 INCL_DEBUG("PiN interaction: one Pion channel chosen" << '\n');
760 weight = counterweight;
761 return new PiNToMultiPionsChannel(2,particle1, particle2);
762 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
763 isElastic = false;
764// PiN -> PiN2Pi channel is chosen
765 INCL_DEBUG("PiN interaction: two Pions channel chosen" << '\n');
766 weight = counterweight;
767 return new PiNToMultiPionsChannel(3,particle1, particle2);
768 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
769 isElastic = false;
770// PiN -> PiN3Pi channel is chosen
771 INCL_DEBUG("PiN interaction: three Pions channel chosen" << '\n');
772 weight = counterweight;
773 return new PiNToMultiPionsChannel(4,particle1, particle2);
774 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX > rChannel) {
775 isElastic = false;
776// PiN -> EtaN channel is chosen
777 INCL_DEBUG("PiN interaction: Eta channel chosen" << '\n');
778 weight = counterweight;
779 return new PiNToEtaChannel(particle1, particle2);
780 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX > rChannel) {
781 isElastic = false;
782// PiN -> OmegaN channel is chosen
783 INCL_DEBUG("PiN interaction: Omega channel chosen" << '\n');
784 weight = counterweight;
785 return new PiNToOmegaChannel(particle1, particle2);
786 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
787 + LKProdCX > rChannel) {
788 isElastic = false;
789 isStrangeProduction = true;
790// PiN -> LK channel is chosen
791 INCL_DEBUG("PiN interaction: LK channel chosen" << '\n');
792 weight = bias_apply;
793 return new NpiToLKChannel(particle1, particle2);
794 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
795 + LKProdCX + SKProdCX > rChannel) {
796 isElastic = false;
797 isStrangeProduction = true;
798// PiN -> SK channel is chosen
799 INCL_DEBUG("PiN interaction: SK channel chosen" << '\n');
800 weight = bias_apply;
801 return new NpiToSKChannel(particle1, particle2);
802 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
803 + LKProdCX + SKProdCX + LKpiProdCX > rChannel) {
804 isElastic = false;
805 isStrangeProduction = true;
806// PiN -> LKpi channel is chosen
807 INCL_DEBUG("PiN interaction: LKpi channel chosen" << '\n');
808 weight = bias_apply;
809 return new NpiToLKpiChannel(particle1, particle2);
810 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
811 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX > rChannel) {
812 isElastic = false;
813 isStrangeProduction = true;
814// PiN -> SKpi channel is chosen
815 INCL_DEBUG("PiN interaction: SKpi channel chosen" << '\n');
816 weight = bias_apply;
817 return new NpiToSKpiChannel(particle1, particle2);
818 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
819 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX > rChannel) {
820 isElastic = false;
821 isStrangeProduction = true;
822// PiN -> LK2pi channel is chosen
823 INCL_DEBUG("PiN interaction: LK2pi channel chosen" << '\n');
824 weight = bias_apply;
825 return new NpiToLK2piChannel(particle1, particle2);
826 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
827 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX > rChannel) {
828 isElastic = false;
829 isStrangeProduction = true;
830// PiN -> SK2pi channel is chosen
831 INCL_DEBUG("PiN interaction: SK2pi channel chosen" << '\n');
832 weight = bias_apply;
833 return new NpiToSK2piChannel(particle1, particle2);
834 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
835 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX > rChannel) {
836 isElastic = false;
837 isStrangeProduction = true;
838// PiN -> NKKb channel is chosen
839 INCL_DEBUG("PiN interaction: NKKb channel chosen" << '\n');
840 weight = bias_apply;
841 return new NpiToNKKbChannel(particle1, particle2);
842 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
843 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX> rChannel) {
844 isElastic = false;
845 isStrangeProduction = true;
846// PiN -> Missinge Strangeness channel is chosen
847 INCL_DEBUG("PiN interaction: Missinge Strangeness channel chosen" << '\n');
848 weight = bias_apply;
849 return new NpiToMissingStrangenessChannel(particle1, particle2);
850 }
851 else {
852 INCL_WARN("inconsistency within the PiN Cross Sections (sum!=inelastic)" << '\n');
853 if(MissingCX>0.) {
854 INCL_WARN("Returning a Missinge Strangeness channel" << '\n');
855 weight = bias_apply;
856 isElastic = false;
857 isStrangeProduction = true;
858 return new NpiToMissingStrangenessChannel(particle1, particle2);
859 } else if(NKKbProdCX>0.) {
860 INCL_WARN("Returning a NKKb channel" << '\n');
861 weight = bias_apply;
862 isElastic = false;
863 isStrangeProduction = true;
864 return new NpiToNKKbChannel(particle1, particle2);
865 } else if(SK2piProdCX>0.) {
866 INCL_WARN("Returning a SK2pi channel" << '\n');
867 weight = bias_apply;
868 isElastic = false;
869 isStrangeProduction = true;
870 return new NpiToSK2piChannel(particle1, particle2);
871 } else if(LK2piProdCX>0.) {
872 INCL_WARN("Returning a LK2pi channel" << '\n');
873 weight = bias_apply;
874 isElastic = false;
875 isStrangeProduction = true;
876 return new NpiToLK2piChannel(particle1, particle2);
877 } else if(SKpiProdCX>0.) {
878 INCL_WARN("Returning a SKpi channel" << '\n');
879 weight = bias_apply;
880 isElastic = false;
881 isStrangeProduction = true;
882 return new NpiToSKpiChannel(particle1, particle2);
883 } else if(LKpiProdCX>0.) {
884 INCL_WARN("Returning a LKpi channel" << '\n');
885 weight = bias_apply;
886 isElastic = false;
887 isStrangeProduction = true;
888 return new NpiToLKpiChannel(particle1, particle2);
889 } else if(SKProdCX>0.) {
890 INCL_WARN("Returning a SK channel" << '\n');
891 weight = bias_apply;
892 isElastic = false;
893 isStrangeProduction = true;
894 return new NpiToSKChannel(particle1, particle2);
895 } else if(LKProdCX>0.) {
896 INCL_WARN("Returning a LK channel" << '\n');
897 weight = bias_apply;
898 isElastic = false;
899 isStrangeProduction = true;
900 return new NpiToLKChannel(particle1, particle2);
901 } else if(omegaProductionCX>0.) {
902 INCL_WARN("Returning a Omega channel" << '\n');
903 weight = counterweight;
904 isElastic = false;
905 return new PiNToOmegaChannel(particle1, particle2);
906 } else if(etaProductionCX>0.) {
907 INCL_WARN("Returning a Eta channel" << '\n');
908 weight = counterweight;
909 isElastic = false;
910 return new PiNToEtaChannel(particle1, particle2);
911 } else if(threePiProductionCX>0.) {
912 INCL_WARN("Returning a 3pi channel" << '\n');
913 weight = counterweight;
914 isElastic = false;
915 return new PiNToMultiPionsChannel(4,particle1, particle2);
916 } else if(twoPiProductionCX>0.) {
917 INCL_WARN("Returning a 2pi channel" << '\n');
918 weight = counterweight;
919 isElastic = false;
920 return new PiNToMultiPionsChannel(3,particle1, particle2);
921 } else if(onePiProductionCX>0.) {
922 INCL_WARN("Returning a 1pi channel" << '\n');
923 weight = counterweight;
924 isElastic = false;
925 return new PiNToMultiPionsChannel(2,particle1, particle2);
926 } else if(deltaProductionCX>0.) {
927 INCL_WARN("Returning a delta-production channel" << '\n');
928 weight = counterweight;
929 isElastic = false;
930 return new PiNToDeltaChannel(particle1, particle2);
931 } else {
932 INCL_WARN("Returning an elastic channel" << '\n');
933 weight = counterweight;
934 isElastic = true;
935 return new PiNElasticChannel(particle1, particle2);
936 }
937 }
938 } else if ((particle1->isNucleon() && particle2->isEta()) || (particle2->isNucleon() && particle1->isEta())) {
939//// EtaN
940
941 const G4double elasticCX = CrossSections::elastic(particle1, particle2);
942 const G4double onePiProductionCX = CrossSections::etaNToPiN(particle1, particle2);
943 const G4double twoPiProductionCX = CrossSections::etaNToPiPiN(particle1, particle2);
944 const G4double totCX=CrossSections::total(particle1, particle2);
945// assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
946
947 const G4double rChannel=Random::shoot() * totCX;
948
949 if(elasticCX > rChannel) {
950// Elastic EtaN channel
951 isElastic = true;
952 INCL_DEBUG("EtaN interaction: elastic channel chosen" << '\n');
953 return new EtaNElasticChannel(particle1, particle2);
954 } else if(elasticCX + onePiProductionCX > rChannel) {
955 isElastic = false;
956// EtaN -> EtaPiN channel is chosen
957 INCL_DEBUG("EtaN interaction: PiN channel chosen" << '\n');
958 return new EtaNToPiNChannel(particle1, particle2);
959 } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
960 isElastic = false;
961// EtaN -> EtaPiPiN channel is chosen
962 INCL_DEBUG("EtaN interaction: PiPiN channel chosen" << '\n');
963 return new EtaNToPiPiNChannel(particle1, particle2);
964 }
965
966 else {
967 INCL_WARN("inconsistency within the EtaN Cross Sections (sum!=inelastic)" << '\n');
968 if(twoPiProductionCX>0.) {
969 INCL_WARN("Returning a PiPiN channel" << '\n');
970 isElastic = false;
971 return new EtaNToPiPiNChannel(particle1, particle2);
972 } else if(onePiProductionCX>0.) {
973 INCL_WARN("Returning a PiN channel" << '\n');
974 isElastic = false;
975 return new EtaNToPiNChannel(particle1, particle2);
976 } else {
977 INCL_WARN("Returning an elastic channel" << '\n');
978 isElastic = true;
979 return new EtaNElasticChannel(particle1, particle2);
980 }
981 }
982
983 } else if ((particle1->isNucleon() && particle2->isOmega()) || (particle2->isNucleon() && particle1->isOmega())) {
984//// OmegaN
985
986 const G4double elasticCX = CrossSections::elastic(particle1, particle2);
987 const G4double onePiProductionCX = CrossSections::omegaNToPiN(particle1, particle2);
988 const G4double twoPiProductionCX = CrossSections::omegaNToPiPiN(particle1, particle2);
989 const G4double totCX=CrossSections::total(particle1, particle2);
990// assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
991
992 const G4double rChannel=Random::shoot() * totCX;
993
994 if(elasticCX > rChannel) {
995// Elastic OmegaN channel
996 isElastic = true;
997 INCL_DEBUG("OmegaN interaction: elastic channel chosen" << '\n');
998 return new OmegaNElasticChannel(particle1, particle2);
999 } else if(elasticCX + onePiProductionCX > rChannel) {
1000 isElastic = false;
1001// OmegaN -> PiN channel is chosen
1002 INCL_DEBUG("OmegaN interaction: PiN channel chosen" << '\n');
1003 return new OmegaNToPiNChannel(particle1, particle2);
1004 } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
1005 isElastic = false;
1006// OmegaN -> PiPiN channel is chosen
1007 INCL_DEBUG("OmegaN interaction: PiPiN channel chosen" << '\n');
1008 return new OmegaNToPiPiNChannel(particle1, particle2);
1009 }
1010 else {
1011 INCL_WARN("inconsistency within the OmegaN Cross Sections (sum!=inelastic)" << '\n');
1012 if(twoPiProductionCX>0.) {
1013 INCL_WARN("Returning a PiPiN channel" << '\n');
1014 isElastic = false;
1015 return new OmegaNToPiPiNChannel(particle1, particle2);
1016 } else if(onePiProductionCX>0.) {
1017 INCL_WARN("Returning a PiN channel" << '\n');
1018 isElastic = false;
1019 return new OmegaNToPiNChannel(particle1, particle2);
1020 } else {
1021 INCL_WARN("Returning an elastic channel" << '\n');
1022 isElastic = true;
1023 return new OmegaNElasticChannel(particle1, particle2);
1024 }
1025 }
1026 } else if ((particle1->isNucleon() && particle2->isKaon()) || (particle2->isNucleon() && particle1->isKaon())) {
1027//// KN
1028 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1029 const G4double quasielasticCX = CrossSections::NKToNK(particle1,particle2);
1030 const G4double NKToNKpiCX = CrossSections::NKToNKpi(particle1,particle2);
1031 const G4double NKToNK2piCX = CrossSections::NKToNK2pi(particle1,particle2);
1032 const G4double totCX=CrossSections::total(particle1, particle2);
1033// assert(std::fabs(totCX-elasticCX-quasielasticCX-NKToNKpiCX-NKToNK2piCX)<0.1);
1034
1035 const G4double rChannel=Random::shoot() * totCX;
1036 if(elasticCX > rChannel){
1037// Elastic KN channel is chosen
1038 isElastic = true;
1039 INCL_DEBUG("KN interaction: elastic channel chosen" << '\n');
1040 return new NKElasticChannel(particle1, particle2);
1041 } else if(elasticCX + quasielasticCX > rChannel){
1042// Quasi-elastic KN channel is chosen
1043 isElastic = false; // true ??
1044 INCL_DEBUG("KN interaction: quasi-elastic channel chosen" << '\n');
1045 return new NKToNKChannel(particle1, particle2);
1046 } else if(elasticCX + quasielasticCX + NKToNKpiCX > rChannel){
1047// KN -> NKpi channel is chosen
1048 isElastic = false;
1049 INCL_DEBUG("KN interaction: NKpi channel chosen" << '\n');
1050 return new NKToNKpiChannel(particle1, particle2);
1051 } else if(elasticCX + quasielasticCX + NKToNKpiCX + NKToNK2piCX > rChannel){
1052// KN -> NK2pi channel is chosen
1053 isElastic = false;
1054 INCL_DEBUG("KN interaction: NK2pi channel chosen" << '\n');
1055 return new NKToNK2piChannel(particle1, particle2);
1056 } else {
1057 INCL_WARN("inconsistency within the KN Cross Sections (sum!=inelastic)" << '\n');
1058 if(NKToNK2piCX>0.) {
1059 INCL_WARN("Returning a NKToNK2pi channel" << '\n');
1060 isElastic = false;
1061 return new NKToNK2piChannel(particle1, particle2);
1062 } else if(NKToNKpiCX>0.) {
1063 INCL_WARN("Returning a NKToNKpi channel" << '\n');
1064 isElastic = false;
1065 return new NKToNKpiChannel(particle1, particle2);
1066 } else if(quasielasticCX>0.) {
1067 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1068 isElastic = false; // true ??
1069 return new NKToNKChannel(particle1, particle2);
1070 } else {
1071 INCL_WARN("Returning an elastic channel" << '\n');
1072 isElastic = true;
1073 return new NKElasticChannel(particle1, particle2);
1074 }
1075 }
1076 } else if ((particle1->isNucleon() && particle2->isAntiKaon()) || (particle2->isNucleon() && particle1->isAntiKaon())) {
1077//// KbN
1078 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1079 const G4double quasielasticCX = CrossSections::NKbToNKb(particle1,particle2);
1080 const G4double NKbToNKbpiCX = CrossSections::NKbToNKbpi(particle1,particle2);
1081 const G4double NKbToNKb2piCX = CrossSections::NKbToNKb2pi(particle1,particle2);
1082 const G4double NKbToLpiCX = CrossSections::NKbToLpi(particle1,particle2);
1083 const G4double NKbToL2piCX = CrossSections::NKbToL2pi(particle1,particle2);
1084 const G4double NKbToSpiCX = CrossSections::NKbToSpi(particle1,particle2);
1085 const G4double NKbToS2piCX = CrossSections::NKbToS2pi(particle1,particle2);
1086 const G4double totCX=CrossSections::total(particle1, particle2);
1087// assert(std::fabs(totCX-elasticCX-quasielasticCX-NKbToNKbpiCX-NKbToNKb2piCX-NKbToLpiCX-NKbToL2piCX-NKbToSpiCX-NKbToS2piCX)<0.1);
1088
1089 const G4double rChannel=Random::shoot() * totCX;
1090 if(elasticCX > rChannel){
1091// Elastic KbN channel is chosen
1092 isElastic = true;
1093 INCL_DEBUG("KbN interaction: elastic channel chosen" << '\n');
1094 return new NKbElasticChannel(particle1, particle2);
1095 } else if(elasticCX + quasielasticCX > rChannel){
1096// Quasi-elastic KbN channel is chosen
1097 isElastic = false; // true ??
1098 INCL_DEBUG("KbN interaction: quasi-elastic channel chosen" << '\n');
1099 return new NKbToNKbChannel(particle1, particle2);
1100 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX > rChannel){
1101// KbN -> NKbpi channel is chosen
1102 isElastic = false;
1103 INCL_DEBUG("KbN interaction: NKbpi channel chosen" << '\n');
1104 return new NKbToNKbpiChannel(particle1, particle2);
1105 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX > rChannel){
1106// KbN -> NKb2pi channel is chosen
1107 isElastic = false;
1108 INCL_DEBUG("KbN interaction: NKb2pi channel chosen" << '\n');
1109 return new NKbToNKb2piChannel(particle1, particle2);
1110 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX > rChannel){
1111// KbN -> Lpi channel is chosen
1112 isElastic = false;
1113 INCL_DEBUG("KbN interaction: Lpi channel chosen" << '\n');
1114 return new NKbToLpiChannel(particle1, particle2);
1115 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX > rChannel){
1116// KbN -> L2pi channel is chosen
1117 isElastic = false;
1118 INCL_DEBUG("KbN interaction: L2pi channel chosen" << '\n');
1119 return new NKbToL2piChannel(particle1, particle2);
1120 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX > rChannel){
1121// KbN -> Spi channel is chosen
1122 isElastic = false;
1123 INCL_DEBUG("KbN interaction: Spi channel chosen" << '\n');
1124 return new NKbToSpiChannel(particle1, particle2);
1125 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX + NKbToS2piCX > rChannel){
1126// KbN -> S2pi channel is chosen
1127 isElastic = false;
1128 INCL_DEBUG("KbN interaction: S2pi channel chosen" << '\n');
1129 return new NKbToS2piChannel(particle1, particle2);
1130 } else {
1131 INCL_WARN("inconsistency within the KbN Cross Sections (sum!=inelastic)" << '\n');
1132 if(NKbToS2piCX>0.) {
1133 INCL_WARN("Returning a NKbToS2pi channel" << '\n');
1134 isElastic = false;
1135 return new NKbToS2piChannel(particle1, particle2);
1136 } else if(NKbToSpiCX>0.) {
1137 INCL_WARN("Returning a NKbToSpi channel" << '\n');
1138 isElastic = false;
1139 return new NKbToSpiChannel(particle1, particle2);
1140 } else if(NKbToL2piCX>0.) {
1141 INCL_WARN("Returning a NKbToL2pi channel" << '\n');
1142 isElastic = false;
1143 return new NKbToL2piChannel(particle1, particle2);
1144 } else if(NKbToLpiCX>0.) {
1145 INCL_WARN("Returning a NKbToLpi channel" << '\n');
1146 isElastic = false;
1147 return new NKbToLpiChannel(particle1, particle2);
1148 } else if(NKbToNKb2piCX>0.) {
1149 INCL_WARN("Returning a NKbToNKb2pi channel" << '\n');
1150 isElastic = false;
1151 return new NKbToNKb2piChannel(particle1, particle2);
1152 } else if(NKbToNKbpiCX>0.) {
1153 INCL_WARN("Returning a NKbToNKbpi channel" << '\n');
1154 isElastic = false;
1155 return new NKbToNKbpiChannel(particle1, particle2);
1156 } else if(quasielasticCX>0.) {
1157 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1158 isElastic = false; // true ??
1159 return new NKbToNKbChannel(particle1, particle2);
1160 } else {
1161 INCL_WARN("Returning an elastic channel" << '\n');
1162 isElastic = true;
1163 return new NKbElasticChannel(particle1, particle2);
1164 }
1165 }
1166 } else if ((particle1->isNucleon() && particle2->isLambda()) || (particle2->isNucleon() && particle1->isLambda())) {
1167//// NLambda
1168 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1169 const G4double NLToNSCX = CrossSections::NLToNS(particle1,particle2);
1170 const G4double totCX=CrossSections::total(particle1, particle2);
1171// assert(std::fabs(totCX-elasticCX-NLToNSCX)<0.1);
1172
1173 const G4double rChannel=Random::shoot() * totCX;
1174 if(elasticCX > rChannel){
1175// Elastic NLambda channel is chosen
1176 isElastic = true;
1177 INCL_DEBUG("NLambda interaction: elastic channel chosen" << '\n');
1178 return new NYElasticChannel(particle1, particle2);
1179 } else if(elasticCX + NLToNSCX > rChannel){
1180// Quasi-elastic NLambda channel is chosen
1181 isElastic = false; // true ??
1182 INCL_DEBUG("NLambda interaction: quasi-elastic channel chosen" << '\n');
1183 return new NLToNSChannel(particle1, particle2);
1184 } else {
1185 INCL_WARN("inconsistency within the NLambda Cross Sections (sum!=inelastic)" << '\n');
1186 if(NLToNSCX>0.) {
1187 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1188 isElastic = false; // true ??
1189 return new NLToNSChannel(particle1, particle2);
1190 } else {
1191 INCL_WARN("Returning an elastic channel" << '\n');
1192 isElastic = true;
1193 return new NYElasticChannel(particle1, particle2);
1194 }
1195 }
1196 } else if ((particle1->isNucleon() && particle2->isSigma()) || (particle2->isNucleon() && particle1->isSigma())) {
1197//// NSigma
1198 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1199 const G4double NSToNLCX = CrossSections::NSToNL(particle1,particle2);
1200 const G4double NSToNSCX = CrossSections::NSToNS(particle1,particle2);
1201 const G4double totCX=CrossSections::total(particle1, particle2);
1202// assert(std::fabs(totCX-elasticCX-NSToNLCX-NSToNSCX)<0.1);
1203
1204 const G4double rChannel=Random::shoot() * totCX;
1205 if(elasticCX > rChannel){
1206// Elastic NSigma channel is chosen
1207 isElastic = true;
1208 INCL_DEBUG("NSigma interaction: elastic channel chosen" << '\n');
1209 return new NYElasticChannel(particle1, particle2);
1210 } else if(elasticCX + NSToNLCX > rChannel){
1211// NSigma -> NLambda channel is chosen
1212 isElastic = false; // true ??
1213 INCL_DEBUG("NSigma interaction: NLambda channel chosen" << '\n');
1214 return new NSToNLChannel(particle1, particle2);
1215 } else if(elasticCX + NSToNLCX + NSToNSCX > rChannel){
1216// NSigma -> NSigma quasi-elastic channel is chosen
1217 isElastic = false; // true ??
1218 INCL_DEBUG("NSigma interaction: NSigma quasi-elastic channel chosen" << '\n');
1219 return new NSToNSChannel(particle1, particle2);
1220 } else {
1221 INCL_WARN("inconsistency within the NSigma Cross Sections (sum!=inelastic)" << '\n');
1222 if(NSToNSCX>0.) {
1223 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1224 isElastic = false; // true ??
1225 return new NSToNSChannel(particle1, particle2);
1226 } else if(NSToNLCX>0.) {
1227 INCL_WARN("Returning a NLambda channel" << '\n');
1228 isElastic = false; // true ??
1229 return new NSToNLChannel(particle1, particle2);
1230 } else {
1231 INCL_WARN("Returning an elastic channel" << '\n');
1232 isElastic = true;
1233 return new NYElasticChannel(particle1, particle2);
1234 }
1235 }
1236 }
1237
1238 else {
1239 INCL_DEBUG("BinaryCollisionAvatar can only handle nucleons (for the moment)."
1240 << '\n'
1241 << particle1->print()
1242 << '\n'
1243 << particle2->print()
1244 << '\n');
1245 InteractionAvatar::restoreParticles();
1246 return NULL;
1247 }
1248 }
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
G4INCL::Book::getAcceptedCollisions
G4int getAcceptedCollisions() const
Definition: G4INCLBook.hh:100
G4INCL::InteractionAvatar::restoreParticles
void restoreParticles() const
Restore the state of both particles.
Definition: G4INCLInteractionAvatar.cc:340
G4INCL::Nucleus::getStore
Store * getStore() const
Definition: G4INCLNucleus.hh:347
G4INCL::Particle::isLambda
G4bool isLambda() const
Is this a Lambda?
Definition: G4INCLParticle.hh:379
G4INCL::Particle::isOmega
G4bool isOmega() const
Is this an omega?
Definition: G4INCLParticle.hh:353
G4INCL::Particle::isSigma
G4bool isSigma() const
Is this a Sigma?
Definition: G4INCLParticle.hh:370
G4INCL::Particle::getPosition
const G4INCL::ThreeVector & getPosition() const
Definition: G4INCLParticle.hh:820
G4INCL::Particle::MergeVectorBias
static std::vector< G4int > MergeVectorBias(Particle const *const p1, Particle const *const p2)
Definition: G4INCLParticle.cc:229
G4INCL::Particle::isEta
G4bool isEta() const
Is this an eta?
Definition: G4INCLParticle.hh:350
G4INCL::Particle::isAntiKaon
G4bool isAntiKaon() const
Is this an antiKaon?
Definition: G4INCLParticle.hh:376
G4INCL::Particle::isKaon
G4bool isKaon() const
Is this a Kaon?
Definition: G4INCLParticle.hh:373
G4INCL::Particle::print
std::string print() const
Definition: G4INCLParticle.hh:957
G4INCL::Particle::isDelta
G4bool isDelta() const
Is it a Delta?
Definition: G4INCLParticle.hh:365
G4INCL::Particle::getBiasFromVector
static G4double getBiasFromVector(std::vector< G4int > VectorBias)
Definition: G4INCLParticle.cc:217
G4INCL::Particle::isNucleon
G4bool isNucleon() const
Definition: G4INCLParticle.hh:307
G4INCL::Store::getBook
Book & getBook()
Definition: G4INCLStore.hh:259
G4INCL::ThreeVector::dot
G4double dot(const ThreeVector &v) const
Definition: G4INCLThreeVector.hh:97
G4INCL::ThreeVector::mag2
G4double mag2() const
Definition: G4INCLThreeVector.hh:78
G4INCL::CrossSections::NSToNL
G4double NSToNL(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:259
G4INCL::CrossSections::NNToNNKKb
G4double NNToNNKKb(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:186
G4INCL::CrossSections::elastic
G4double elastic(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:56
G4INCL::CrossSections::NpiToSK
G4double NpiToSK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:215
G4INCL::CrossSections::piNToOmegaN
G4double piNToOmegaN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:88
G4INCL::CrossSections::NKbToNKb2pi
G4double NKbToNKb2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:303
G4INCL::CrossSections::NNToNDeltaOmega
G4double NNToNDeltaOmega(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:145
G4INCL::CrossSections::NDeltaToDeltaLK
G4double NDeltaToDeltaLK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:200
G4INCL::CrossSections::NNToNSKpi
G4double NNToNSKpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:174
G4INCL::CrossSections::etaNToPiN
G4double etaNToPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:96
G4INCL::CrossSections::NDeltaToNNKKb
G4double NDeltaToNNKKb(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:207
G4INCL::CrossSections::NKbToLpi
G4double NKbToLpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:287
G4INCL::CrossSections::NNToNLK2pi
G4double NNToNLK2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:178
G4INCL::CrossSections::etaNToPiPiN
G4double etaNToPiPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:100
G4INCL::CrossSections::NKbToS2pi
G4double NKbToS2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:291
G4INCL::CrossSections::piNToEtaN
G4double piNToEtaN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:84
G4INCL::CrossSections::omegaNToPiN
G4double omegaNToPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:104
G4INCL::CrossSections::NDeltaToNSK
G4double NDeltaToNSK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:197
G4INCL::CrossSections::NKbToSpi
G4double NKbToSpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:283
G4INCL::CrossSections::NNToNDelta
G4double NNToNDelta(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:68
G4INCL::CrossSections::NNToNSK
G4double NNToNSK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:166
G4INCL::CrossSections::NDeltaToNLK
G4double NDeltaToNLK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:194
G4INCL::CrossSections::NLToNS
G4double NLToNS(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:255
G4INCL::CrossSections::piNToDelta
G4double piNToDelta(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:76
G4INCL::CrossSections::NKbToNKb
G4double NKbToNKb(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:279
G4INCL::CrossSections::NNToNSK2pi
G4double NNToNSK2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:182
G4INCL::CrossSections::NNToNLK
G4double NNToNLK(Particle const *const p1, Particle const *const p2)
Strange cross sections.
Definition: G4INCLCrossSections.cc:162
G4INCL::CrossSections::NNToNNEtaExclu
G4double NNToNNEtaExclu(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:120
G4INCL::CrossSections::NNToNDeltaEta
G4double NNToNDeltaEta(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:128
G4INCL::CrossSections::NNToNNOmegaxPi
G4double NNToNNOmegaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:141
G4INCL::CrossSections::NNToNNEtaxPi
G4double NNToNNEtaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:124
G4INCL::CrossSections::NSToNS
G4double NSToNS(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:263
G4INCL::CrossSections::NpiToNKKb
G4double NpiToNKKb(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:247
G4INCL::CrossSections::NpiToLK2pi
G4double NpiToLK2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:239
G4INCL::CrossSections::omegaNToPiPiN
G4double omegaNToPiPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:108
G4INCL::CrossSections::NDeltaToDeltaSK
G4double NDeltaToDeltaSK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:203
G4INCL::CrossSections::NpiToMissingStrangeness
G4double NpiToMissingStrangeness(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:251
G4INCL::CrossSections::total
G4double total(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:60
G4INCL::CrossSections::NDeltaToNN
G4double NDeltaToNN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:64
G4INCL::CrossSections::NpiToLKpi
G4double NpiToLKpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:231
G4INCL::CrossSections::NKbToNKbpi
G4double NKbToNKbpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:299
G4INCL::CrossSections::NNToNLKpi
G4double NNToNLKpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:170
G4INCL::CrossSections::NKbToL2pi
G4double NKbToL2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:295
G4INCL::CrossSections::NNToxPiNN
G4double NNToxPiNN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:72
G4INCL::CrossSections::NKToNK2pi
G4double NKToNK2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:275
G4INCL::CrossSections::NKToNKpi
G4double NKToNKpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:271
G4INCL::CrossSections::NNToNNOmegaExclu
G4double NNToNNOmegaExclu(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:137
G4INCL::CrossSections::NpiToSK2pi
G4double NpiToSK2pi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:243
G4INCL::CrossSections::NKToNK
G4double NKToNK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:267
G4INCL::CrossSections::NNToMissingStrangeness
G4double NNToMissingStrangeness(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:190
G4INCL::CrossSections::piNToxPiN
G4double piNToxPiN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:80
G4INCL::CrossSections::NpiToLK
G4double NpiToLK(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:211
G4INCL::CrossSections::NpiToSKpi
G4double NpiToSKpi(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:235
G4INCL::KinematicsUtils::squareTotalEnergyInCM
G4double squareTotalEnergyInCM(Particle const *const p1, Particle const *const p2)
Definition: G4INCLKinematicsUtils.cc:98
G4INCL::Math::tenPi
const G4double tenPi
Definition: G4INCLGlobals.hh:70
G4INCL::Random::shoot
G4double shoot()
Definition: G4INCLRandom.cc:93

◆ getCutNN()

static G4double G4INCL::BinaryCollisionAvatar::getCutNN ( )
inlinestatic

Definition at line 78 of file G4INCLBinaryCollisionAvatar.hh.

78{ return cutNN; }

◆ getCutNNSquared()

static G4double G4INCL::BinaryCollisionAvatar::getCutNNSquared ( )
inlinestatic

Definition at line 80 of file G4INCLBinaryCollisionAvatar.hh.

80{ return cutNNSquared; }

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

◆ getParticles()

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

Implements G4INCL::IAvatar.

Definition at line 61 of file G4INCLBinaryCollisionAvatar.hh.

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

◆ postInteraction()

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

Implements G4INCL::IAvatar.

Definition at line 1256 of file G4INCLBinaryCollisionAvatar.cc.

1256 {
1257 // Call the postInteraction method of the parent class
1258 // (provides Pauli blocking and enforces energy conservation)
1259 InteractionAvatar::postInteraction(fs);
1260
1261 switch(fs->getValidity()) {
1262 case PauliBlockedFS:
1263 theNucleus->getStore()->getBook().incrementBlockedCollisions();
1264 break;
1265 case NoEnergyConservationFS:
1266 case ParticleBelowFermiFS:
1267 case ParticleBelowZeroFS:
1268 break;
1269 case ValidFS:
1270 Book &theBook = theNucleus->getStore()->getBook();
1271 theBook.incrementAcceptedCollisions();
1272
1273 if(theBook.getAcceptedCollisions() == 1) {
1274 // Store time and cross section of the first collision
1275 G4double t = theBook.getCurrentTime();
1276 theBook.setFirstCollisionTime(t);
1277 theBook.setFirstCollisionXSec(oldXSec);
1278 // Increase the number of Kaon by 1
1279 if(isStrangeProduction) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()+1);
1280 // Store position and momentum of the spectator on the first
1281 // collision
1282 if((isParticle1Spectator && isParticle2Spectator) || (!isParticle1Spectator && !isParticle2Spectator)) {
1283 INCL_ERROR("First collision must be within a target spectator and a non-target spectator");
1284 }
1285 if(isParticle1Spectator) {
1286 theBook.setFirstCollisionSpectatorPosition(backupParticle1->getPosition().mag());
1287 theBook.setFirstCollisionSpectatorMomentum(backupParticle1->getMomentum().mag());
1288 } else {
1289 theBook.setFirstCollisionSpectatorPosition(backupParticle2->getPosition().mag());
1290 theBook.setFirstCollisionSpectatorMomentum(backupParticle2->getMomentum().mag());
1291 }
1292
1293 // Store the elasticity of the first collision
1294 theBook.setFirstCollisionIsElastic(isElastic);
1295 }
1296 }
1297 return;
1298 }
G4INCL::Book::incrementAcceptedCollisions
void incrementAcceptedCollisions()
Definition: G4INCLBook.hh:72
G4INCL::Book::incrementBlockedCollisions
void incrementBlockedCollisions()
Definition: G4INCLBook.hh:73
G4INCL::InteractionAvatar::backupParticle2
static G4ThreadLocal Particle * backupParticle2
Definition: G4INCLInteractionAvatar.hh:108
G4INCL::InteractionAvatar::backupParticle1
static G4ThreadLocal Particle * backupParticle1
Definition: G4INCLInteractionAvatar.hh:108
G4INCL::Particle::setNumberOfKaon
void setNumberOfKaon(const G4int NK)
Definition: G4INCLParticle.hh:1058
G4INCL::Particle::getMomentum
const G4INCL::ThreeVector & getMomentum() const
Definition: G4INCLParticle.hh:798
G4INCL::Particle::getNumberOfKaon
G4int getNumberOfKaon() const
Number of Kaon inside de nucleus.
Definition: G4INCLParticle.hh:1057
G4INCL::ThreeVector::mag
G4double mag() const
Definition: G4INCLThreeVector.hh:73
G4INCL::ParticleBelowZeroFS
@ ParticleBelowZeroFS
Definition: G4INCLFinalState.hh:58
G4INCL::NoEnergyConservationFS
@ NoEnergyConservationFS
Definition: G4INCLFinalState.hh:56
G4INCL::ParticleBelowFermiFS
@ ParticleBelowFermiFS
Definition: G4INCLFinalState.hh:57

◆ preInteraction()

void G4INCL::BinaryCollisionAvatar::preInteraction ( )
virtual

Implements G4INCL::IAvatar.

Definition at line 1250 of file G4INCLBinaryCollisionAvatar.cc.

1250 {
1251 isParticle1Spectator = particle1->isTargetSpectator();
1252 isParticle2Spectator = particle2->isTargetSpectator();
1253 InteractionAvatar::preInteraction();
1254 }
G4INCL::Particle::isTargetSpectator
G4bool isTargetSpectator() const
Definition: G4INCLParticle.hh:326

◆ setBias()

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

Set the global bias factor.

Definition at line 86 of file G4INCLBinaryCollisionAvatar.hh.

86{ bias=b; }

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

◆ setCutNN()

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

Definition at line 73 of file G4INCLBinaryCollisionAvatar.hh.

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

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

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