Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCLBinaryCollisionAvatar.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// INCL++ intra-nuclear cascade model
27// Alain Boudard, CEA-Saclay, France
28// Joseph Cugnon, University of Liege, Belgium
29// Jean-Christophe David, CEA-Saclay, France
30// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31// Sylvie Leray, CEA-Saclay, France
32// Davide Mancusi, CEA-Saclay, France
33//
34#define INCLXX_IN_GEANT4_MODE 1
35
36#include "globals.hh"
37
38/*
39 * G4INCLBinaryCollisionAvatar.cc
40 *
41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi
43 */
44
45#include "G4INCLBinaryCollisionAvatar.hh"
46#include "G4INCLElasticChannel.hh"
47#include "G4INCLRecombinationChannel.hh"
48#include "G4INCLDeltaProductionChannel.hh"
49#include "G4INCLNNToMultiPionsChannel.hh"
50#include "G4INCLNNToNNEtaChannel.hh"
51#include "G4INCLNDeltaEtaProductionChannel.hh"
52#include "G4INCLNNEtaToMultiPionsChannel.hh"
53#include "G4INCLNNToNNOmegaChannel.hh"
54#include "G4INCLNDeltaOmegaProductionChannel.hh"
55#include "G4INCLNNOmegaToMultiPionsChannel.hh"
56#include "G4INCLCrossSections.hh"
57#include "G4INCLKinematicsUtils.hh"
58#include "G4INCLRandom.hh"
59#include "G4INCLParticleTable.hh"
60#include "G4INCLPauliBlocking.hh"
61#include "G4INCLPiNElasticChannel.hh"
62#include "G4INCLPiNToDeltaChannel.hh"
63#include "G4INCLPiNToMultiPionsChannel.hh"
64#include "G4INCLPiNToEtaChannel.hh"
65#include "G4INCLPiNToOmegaChannel.hh"
66#include "G4INCLEtaNElasticChannel.hh"
67#include "G4INCLEtaNToPiNChannel.hh"
68#include "G4INCLEtaNToPiPiNChannel.hh"
69#include "G4INCLOmegaNElasticChannel.hh"
70#include "G4INCLOmegaNToPiNChannel.hh"
71#include "G4INCLEtaOrOmegaNToLKChannel.hh"
72#include "G4INCLEtaOrOmegaNToSKChannel.hh"
73#include "G4INCLNNToNLKChannel.hh"
74#include "G4INCLNNToNSKChannel.hh"
75#include "G4INCLNNToNLKpiChannel.hh"
76#include "G4INCLNNToNSKpiChannel.hh"
77#include "G4INCLNNToNLK2piChannel.hh"
78#include "G4INCLNNToNSK2piChannel.hh"
79#include "G4INCLNNToNNKKbChannel.hh"
80#include "G4INCLNNToMissingStrangenessChannel.hh"
81#include "G4INCLNDeltaToNLKChannel.hh"
82#include "G4INCLNDeltaToNSKChannel.hh"
83#include "G4INCLNDeltaToDeltaLKChannel.hh"
84#include "G4INCLNDeltaToDeltaSKChannel.hh"
85#include "G4INCLNDeltaToNNKKbChannel.hh"
86#include "G4INCLNpiToLKChannel.hh"
87#include "G4INCLNpiToSKChannel.hh"
88#include "G4INCLNpiToLKpiChannel.hh"
89#include "G4INCLNpiToSKpiChannel.hh"
90#include "G4INCLNpiToLK2piChannel.hh"
91#include "G4INCLNpiToSK2piChannel.hh"
92#include "G4INCLNpiToNKKbChannel.hh"
93#include "G4INCLNpiToMissingStrangenessChannel.hh"
94#include "G4INCLNKElasticChannel.hh"
95#include "G4INCLNKToNKChannel.hh"
96#include "G4INCLNKToNKpiChannel.hh"
97#include "G4INCLNKToNK2piChannel.hh"
98#include "G4INCLNKbElasticChannel.hh"
99#include "G4INCLNKbToNKbChannel.hh"
100#include "G4INCLNKbToNKbpiChannel.hh"
101#include "G4INCLNKbToNKb2piChannel.hh"
102#include "G4INCLNKbToLpiChannel.hh"
103#include "G4INCLNKbToL2piChannel.hh"
104#include "G4INCLNKbToSpiChannel.hh"
105#include "G4INCLNKbToS2piChannel.hh"
106#include "G4INCLNYElasticChannel.hh"
107#include "G4INCLNLToNSChannel.hh"
108#include "G4INCLNSToNLChannel.hh"
109#include "G4INCLNSToNSChannel.hh"
110#include "G4INCLOmegaNToPiPiNChannel.hh"
111#include "G4INCLStore.hh"
112#include "G4INCLBook.hh"
113#include "G4INCLLogger.hh"
114#include <string>
115#include <sstream>
116// #include <cassert>
117#include "G4INCLNNbarElasticChannel.hh"
118#include "G4INCLNNbarCEXChannel.hh"
119#include "G4INCLNNbarToLLbarChannel.hh"
120#include "G4INCLNNbarToNNbarpiChannel.hh"
121#include "G4INCLNNbarToNNbar2piChannel.hh"
122#include "G4INCLNNbarToNNbar3piChannel.hh"
123#include "G4INCLNNbarToAnnihilationChannel.hh"
124
125namespace G4INCL {
126
127 // WARNING: if you update the default cutNN value, make sure you update the
128 // cutNNSquared variable, too.
129 G4ThreadLocal G4double BinaryCollisionAvatar::cutNN = 1910.0;
130 G4ThreadLocal G4double BinaryCollisionAvatar::cutNNSquared = 3648100.0; // 1910.0 * 1910.0
131 G4ThreadLocal G4double BinaryCollisionAvatar::bias = 1.;
132
133 BinaryCollisionAvatar::BinaryCollisionAvatar(G4double time, G4double crossSection,
134 G4INCL::Nucleus *n, G4INCL::Particle *p1, G4INCL::Particle *p2)
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 }
143
146
147void BinaryCollisionAvatar::generateSrcPairsMethod(ParticleList &theList,
148 const G4int theN,
149 const G4int theZ,
150 G4INCL::Particle *p1,
151 G4INCL::Particle *p2) {
152
153 std::vector<ThreeVector> posp;
154 std::vector<ThreeVector> posn;
155 posp.resize(theZ);
156 posn.resize(theN);
157
158 for (int i = 0; i < theZ + theN; ++i) {
159 theList[i]->setNumberOfSrcPair(0);
160 }
161 theNucleus->getStore()->getBook().setSrcPairs(0);
162
163 // Check that the loops in theZ and theN do what we wish!
164 G4int npairs = 0;
165 for (G4int i = 0; i < theZ; ++i) {
166 Particle *ap = theList[i];
167 posp[i] = ap->getPosition();
168 for (G4int j = 0; j < theN; ++j) {
169 Particle *an = theList[j + theZ];
170 posn[j] = an->getPosition();
171 if ((posn[j] - posp[i]).mag() < ParticleTable::getsrcPairDistance() &&
172 ap->getSrcPair() == 0 && an->getSrcPair() == 0 &&
173 ap->getType() != an->getType() && ap->isTargetSpectator() == 1 &&
174 an->isTargetSpectator() == 1) {
175 npairs++;
176 theNucleus->getStore()->getBook().incrementSrcPairs();
177 if (an == p1 || an == p2 || ap == p1 || ap == p2) {
178 theList[i]->setNumberOfSrcPair(npairs);
179 theList[j + theZ]->setNumberOfSrcPair(npairs);
180 }
181 }
182 }
183 }
184
185 G4int nbp = 0, nbn = 0;
186 for (G4int i = 0; i < theZ + theN; ++i) {
187
188 if (theList[i]->getSrcPair() > 0 && theList[i]->getType() == Proton)
189 nbp++;
190
191 if (theList[i]->getSrcPair() > 0 && theList[i]->getType() == Neutron)
192 nbn++;
193 }
194
195 if (nbp != nbn) {
196 INCL_DEBUG("Pairs: " << nbp << " " << nbn << '\n');
197 for (G4int i = 0; i < theZ; ++i) {
198 INCL_DEBUG(ParticleTable::getName(theList[i]->getType())
199 << " " << theList[i]->isTargetSpectator() << " "
200 << theList[i]->getSrcPair());
201 }
202 INCL_DEBUG("----------- End ------------------" << '\n');
203 }
204
205 return;
206}
207
208 G4INCL::IChannel* BinaryCollisionAvatar::getChannel() {
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 }
1579
1580 void BinaryCollisionAvatar::preInteraction() {
1581 isParticle1Spectator = particle1->isTargetSpectator();
1582 isParticle2Spectator = particle2->isTargetSpectator();
1583 InteractionAvatar::preInteraction();
1584 }
1585
1586 void BinaryCollisionAvatar::postInteraction(FinalState *fs) {
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 }
1629
1630 std::string BinaryCollisionAvatar::dump() const {
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 }
1639
1640}
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
G4INCLRecombinationChannel.hh
Delta-nucleon recombination channel.
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4INCL::BinaryCollisionAvatar::BinaryCollisionAvatar
BinaryCollisionAvatar(G4double, G4double, G4INCL::Nucleus *, G4INCL::Particle *, G4INCL::Particle *)
Definition G4INCLBinaryCollisionAvatar.cc:133
G4INCL::Book::setFirstCollisionSpectatorMomentum
void setFirstCollisionSpectatorMomentum(const G4double x)
Definition G4INCLBook.hh:95
G4INCL::Book::setFirstCollisionXSec
void setFirstCollisionXSec(const G4double x)
Definition G4INCLBook.hh:89
G4INCL::Book::incrementSrcPairs
void incrementSrcPairs()
Definition G4INCLBook.hh:78
G4INCL::Book::incrementAcceptedCollisions
void incrementAcceptedCollisions()
Definition G4INCLBook.hh:75
G4INCL::Book::setSrcPairs
void setSrcPairs(G4int n)
Definition G4INCLBook.hh:104
G4INCL::Book::setFirstCollisionTime
void setFirstCollisionTime(const G4double t)
Definition G4INCLBook.hh:86
G4INCL::Book::getCurrentTime
G4double getCurrentTime() const
Definition G4INCLBook.hh:102
G4INCL::Book::getAcceptedCollisions
G4int getAcceptedCollisions() const
Definition G4INCLBook.hh:110
G4INCL::Book::setFirstCollisionSpectatorPosition
void setFirstCollisionSpectatorPosition(const G4double x)
Definition G4INCLBook.hh:92
G4INCL::Book::setFirstCollisionIsElastic
void setFirstCollisionIsElastic(const G4bool e)
Definition G4INCLBook.hh:98
G4INCL::Config::getAtrestThreshold
G4double getAtrestThreshold() const
Get the pbar at rest annihilation threshold.
Definition G4INCLConfig.hh:391
G4INCL::FinalState::getValidity
FinalStateValidity getValidity() const
Definition G4INCLFinalState.hh:88
G4INCL::IAvatar::getType
AvatarType getType() const
Definition G4INCLIAvatar.hh:91
G4INCL::IAvatar::setType
void setType(AvatarType t)
Definition G4INCLIAvatar.hh:94
G4INCL::InteractionAvatar::restoreParticles
void restoreParticles() const
Restore the state of both particles.
Definition G4INCLInteractionAvatar.cc:474
G4INCL::InteractionAvatar::InteractionAvatar
InteractionAvatar(G4double, G4INCL::Nucleus *, G4INCL::Particle *)
Definition G4INCLInteractionAvatar.cc:71
G4INCL::InteractionAvatar::backupParticle2
static G4ThreadLocal Particle * backupParticle2
Definition G4INCLInteractionAvatar.hh:123
G4INCL::InteractionAvatar::backupParticle1
static G4ThreadLocal Particle * backupParticle1
Definition G4INCLInteractionAvatar.hh:123
G4INCL::Nucleus::getStore
Store * getStore() const
Definition G4INCLNucleus.hh:428
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::Store::getBook
Book & getBook()
Definition G4INCLStore.hh:277
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: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::getsrcPairDistance
G4float getsrcPairDistance()
Get the distance between src nucleons.
Definition G4INCLParticleTable.cc:1384
G4INCL::ParticleTable::getName
std::string getName(const ParticleType t)
Get the native INCL name of the particle.
Definition G4INCLParticleTable.cc:674
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::CollisionAvatarType
@ CollisionAvatarType
Definition G4INCLIAvatar.hh:67
G4INCL::NoEnergyConservationFS
@ NoEnergyConservationFS
Definition G4INCLFinalState.hh:56
G4INCL::NbarPTypeInFlight
@ NbarPTypeInFlight
Definition G4INCLNucleus.hh:65
G4INCL::NbarNTypeInFlight
@ NbarNTypeInFlight
Definition G4INCLNucleus.hh:65
G4ThreadLocal
#define G4ThreadLocal
Definition tls.hh:77