ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLBinaryCollisionAvatar.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLBinaryCollisionAvatar.cc
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 "G4INCLNNToNLKChannel.hh"
72 #include "G4INCLNNToNSKChannel.hh"
73 #include "G4INCLNNToNLKpiChannel.hh"
74 #include "G4INCLNNToNSKpiChannel.hh"
75 #include "G4INCLNNToNLK2piChannel.hh"
76 #include "G4INCLNNToNSK2piChannel.hh"
77 #include "G4INCLNNToNNKKbChannel.hh"
78 #include "G4INCLNNToMissingStrangenessChannel.hh"
79 #include "G4INCLNDeltaToNLKChannel.hh"
80 #include "G4INCLNDeltaToNSKChannel.hh"
81 #include "G4INCLNDeltaToDeltaLKChannel.hh"
82 #include "G4INCLNDeltaToDeltaSKChannel.hh"
83 #include "G4INCLNDeltaToNNKKbChannel.hh"
84 #include "G4INCLNpiToLKChannel.hh"
85 #include "G4INCLNpiToSKChannel.hh"
86 #include "G4INCLNpiToLKpiChannel.hh"
87 #include "G4INCLNpiToSKpiChannel.hh"
88 #include "G4INCLNpiToLK2piChannel.hh"
89 #include "G4INCLNpiToSK2piChannel.hh"
90 #include "G4INCLNpiToNKKbChannel.hh"
91 #include "G4INCLNpiToMissingStrangenessChannel.hh"
92 #include "G4INCLNKElasticChannel.hh"
93 #include "G4INCLNKToNKChannel.hh"
94 #include "G4INCLNKToNKpiChannel.hh"
95 #include "G4INCLNKToNK2piChannel.hh"
96 #include "G4INCLNKbElasticChannel.hh"
97 #include "G4INCLNKbToNKbChannel.hh"
98 #include "G4INCLNKbToNKbpiChannel.hh"
99 #include "G4INCLNKbToNKb2piChannel.hh"
100 #include "G4INCLNKbToLpiChannel.hh"
101 #include "G4INCLNKbToL2piChannel.hh"
102 #include "G4INCLNKbToSpiChannel.hh"
103 #include "G4INCLNKbToS2piChannel.hh"
104 #include "G4INCLNYElasticChannel.hh"
105 #include "G4INCLNLToNSChannel.hh"
106 #include "G4INCLNSToNLChannel.hh"
107 #include "G4INCLNSToNSChannel.hh"
108 #include "G4INCLOmegaNToPiPiNChannel.hh"
109 #include "G4INCLStore.hh"
110 #include "G4INCLBook.hh"
111 #include "G4INCLLogger.hh"
112 #include <string>
113 #include <sstream>
114 // #include <cassert>
115 
116 namespace G4INCL {
117 
118  // WARNING: if you update the default cutNN value, make sure you update the
119  // cutNNSquared variable, too.
120  G4ThreadLocal G4double BinaryCollisionAvatar::cutNN = 1910.0;
121  G4ThreadLocal G4double BinaryCollisionAvatar::cutNNSquared = 3648100.0; // 1910.0 * 1910.0
122  G4ThreadLocal G4double BinaryCollisionAvatar::bias = 1.;
123 
124  BinaryCollisionAvatar::BinaryCollisionAvatar(G4double time, G4double crossSection,
125  G4INCL::Nucleus *n, G4INCL::Particle *p1, G4INCL::Particle *p2)
126  : InteractionAvatar(time, n, p1, p2), theCrossSection(crossSection),
127  isParticle1Spectator(false),
128  isParticle2Spectator(false),
129  isElastic(false),
130  isStrangeProduction(false)
131  {
132  setType(CollisionAvatarType);
133  }
134 
135  BinaryCollisionAvatar::~BinaryCollisionAvatar() {
136  }
137 
138  G4INCL::IChannel* BinaryCollisionAvatar::getChannel() {
139  // We already check cutNN at avatar creation time, but we have to check it
140  // again here. For composite projectiles, we might have created independent
141  // avatars with no cutNN before any collision took place.
142  if(particle1->isNucleon()
143  && particle2->isNucleon()
144  && theNucleus->getStore()->getBook().getAcceptedCollisions()!=0) {
145  const G4double energyCM2 = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
146  // Below a certain cut value we don't do anything:
147  if(energyCM2 < cutNNSquared) {
148  INCL_DEBUG("CM energy = sqrt(" << energyCM2 << ") MeV < std::sqrt(" << cutNNSquared
149  << ") MeV = cutNN" << "; returning a NULL channel" << '\n');
150  InteractionAvatar::restoreParticles();
151  return NULL;
152  }
153  }
154 
174  ThreeVector minimumDistance = particle1->getPosition();
175  minimumDistance -= particle2->getPosition();
176  const G4double betaDotX = boostVector.dot(minimumDistance);
177  const G4double minDist = Math::tenPi*(minimumDistance.mag2() + betaDotX*betaDotX / (1.-boostVector.mag2()));
178  if(minDist > theCrossSection) {
179  INCL_DEBUG("CM distance of approach is too small: " << minDist << ">" <<
180  theCrossSection <<"; returning a NULL channel" << '\n');
181  InteractionAvatar::restoreParticles();
182  return NULL;
183  }
184 
189  G4double bias_apply = 1.;
190  if(bias != 1.) bias_apply = Particle::getBiasFromVector(Particle::MergeVectorBias(particle1,particle2)) * bias;
191 
193  if(particle1->isNucleon() && particle2->isNucleon()) {
194 
195  G4double NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
196  G4double NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
197  G4double NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
198  G4double NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
199  G4double NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
200  G4double NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
201  G4double NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
202  G4double NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
203 
204  const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NNToNDelta(particle1, particle2) + CrossSections::NNToxPiNN(1,particle1, particle2)
205  + CrossSections::NNToxPiNN(2,particle1, particle2) + CrossSections::NNToxPiNN(3,particle1, particle2) + CrossSections::NNToxPiNN(4,particle1, particle2)
206  + CrossSections::NNToNNEtaExclu(particle1, particle2) + CrossSections::NNToNDeltaEta(particle1, particle2) + CrossSections::NNToNNEtaxPi(1,particle1, particle2)
207  + CrossSections::NNToNNEtaxPi(2,particle1, particle2) + CrossSections::NNToNNEtaxPi(3,particle1, particle2) + CrossSections::NNToNNEtaxPi(4,particle1, particle2)
208  + CrossSections::NNToNNOmegaExclu(particle1, particle2) + CrossSections::NNToNDeltaOmega(particle1, particle2) + CrossSections::NNToNNOmegaxPi(1,particle1, particle2)
209  + CrossSections::NNToNNOmegaxPi(2,particle1, particle2) + CrossSections::NNToNNOmegaxPi(3,particle1, particle2) + CrossSections::NNToNNOmegaxPi(4,particle1, particle2);
210  const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + NLKpiProductionCX + NSKpiProductionCX + NLK2piProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX)/bias_apply;
211 
212  G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
213 
214  if(counterweight < 0.5) {
215  counterweight = 0.5;
216  bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
217  NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
218  NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
219  NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
220  NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
221  NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
222  NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
223  NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
224  NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
225  }
226 
227 
228  const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
229  const G4double deltaProductionCX = CrossSections::NNToNDelta(particle1, particle2)*counterweight;
230  const G4double onePiProductionCX = CrossSections::NNToxPiNN(1,particle1, particle2)*counterweight;
231  const G4double twoPiProductionCX = CrossSections::NNToxPiNN(2,particle1, particle2)*counterweight;
232  const G4double threePiProductionCX = CrossSections::NNToxPiNN(3,particle1, particle2)*counterweight;
233  const G4double fourPiProductionCX = CrossSections::NNToxPiNN(4,particle1, particle2)*counterweight;
234 
235  const G4double etaProductionCX = CrossSections::NNToNNEtaExclu(particle1, particle2)*counterweight;
236  const G4double etadeltaProductionCX = CrossSections::NNToNDeltaEta(particle1, particle2)*counterweight;
237  const G4double etaonePiProductionCX = CrossSections::NNToNNEtaxPi(1,particle1, particle2)*counterweight;
238  const G4double etatwoPiProductionCX = CrossSections::NNToNNEtaxPi(2,particle1, particle2)*counterweight;
239  const G4double etathreePiProductionCX = CrossSections::NNToNNEtaxPi(3,particle1, particle2)*counterweight;
240  const G4double etafourPiProductionCX = CrossSections::NNToNNEtaxPi(4,particle1, particle2)*counterweight;
241  const G4double omegaProductionCX = CrossSections::NNToNNOmegaExclu(particle1, particle2)*counterweight;
242  const G4double omegadeltaProductionCX = CrossSections::NNToNDeltaOmega(particle1, particle2)*counterweight;
243  const G4double omegaonePiProductionCX = CrossSections::NNToNNOmegaxPi(1,particle1, particle2)*counterweight;
244  const G4double omegatwoPiProductionCX = CrossSections::NNToNNOmegaxPi(2,particle1, particle2)*counterweight;
245  const G4double omegathreePiProductionCX = CrossSections::NNToNNOmegaxPi(3,particle1, particle2)*counterweight;
246  const G4double omegafourPiProductionCX = CrossSections::NNToNNOmegaxPi(4,particle1, particle2)*counterweight;
247 
248  const G4double totCX=CrossSections::total(particle1, particle2);
249 
250 // assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-fourPiProductionCX-NLKProductionCX-NSKProductionCX-NLKpiProductionCX-NSKpiProductionCX-NLK2piProductionCX-NSK2piProductionCX-NNKKbProductionCX-NNMissingCX-etaProductionCX-etadeltaProductionCX-etaonePiProductionCX-etatwoPiProductionCX-etathreePiProductionCX-etafourPiProductionCX-omegaProductionCX-omegadeltaProductionCX-omegaonePiProductionCX-omegatwoPiProductionCX-omegathreePiProductionCX-omegafourPiProductionCX) < 0.5);
251 
252  const G4double rChannel=Random::shoot() * totCX;
253 
254  if(elasticCX > rChannel) {
255 // Elastic NN channel
256  isElastic = true;
257  INCL_DEBUG("NN interaction: elastic channel chosen" << '\n');
258  weight = counterweight;
259  return new ElasticChannel(particle1, particle2);
260  } else if((elasticCX + deltaProductionCX) > rChannel) {
261  isElastic = false;
262 // NN -> N Delta channel is chosen
263  INCL_DEBUG("NN interaction: Delta channel chosen" << '\n');
264  weight = counterweight;
265  return new DeltaProductionChannel(particle1, particle2);
266  } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
267  isElastic = false;
268 // NN -> PiNN channel is chosen
269  INCL_DEBUG("NN interaction: one Pion channel chosen" << '\n');
270  weight = counterweight;
271  return new NNToMultiPionsChannel(1,particle1, particle2);
272  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
273  isElastic = false;
274 // NN -> 2PiNN channel is chosen
275  INCL_DEBUG("NN interaction: two Pions channel chosen" << '\n');
276  weight = counterweight;
277  return new NNToMultiPionsChannel(2,particle1, particle2);
278  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
279  isElastic = false;
280 // NN -> 3PiNN channel is chosen
281  INCL_DEBUG("NN interaction: three Pions channel chosen" << '\n');
282  weight = counterweight;
283  return new NNToMultiPionsChannel(3,particle1, particle2);
284  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX > rChannel) {
285  isElastic = false;
286 // NN -> 4PiNN channel is chosen
287  INCL_DEBUG("NN interaction: four Pions channel chosen" << '\n');
288  weight = counterweight;
289  return new NNToMultiPionsChannel(4,particle1, particle2);
290  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
291  + etaProductionCX > rChannel) {
292  isElastic = false;
293 // NN -> NNEta channel is chosen
294  INCL_DEBUG("NN interaction: Eta channel chosen" << '\n');
295  weight = counterweight;
296  return new NNToNNEtaChannel(particle1, particle2);
297  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
298  + etaProductionCX + etadeltaProductionCX > rChannel) {
299  isElastic = false;
300 // NN -> N Delta Eta channel is chosen
301  INCL_DEBUG("NN interaction: Delta Eta channel chosen" << '\n');
302  weight = counterweight;
303  return new NDeltaEtaProductionChannel(particle1, particle2);
304  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
305  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX > rChannel) {
306  isElastic = false;
307 // NN -> EtaPiNN channel is chosen
308  INCL_DEBUG("NN interaction: Eta + one Pion channel chosen" << '\n');
309  weight = counterweight;
310  return new NNEtaToMultiPionsChannel(1,particle1, particle2);
311  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
312  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX > rChannel) {
313  isElastic = false;
314 // NN -> Eta2PiNN channel is chosen
315  INCL_DEBUG("NN interaction: Eta + two Pions channel chosen" << '\n');
316  weight = counterweight;
317  return new NNEtaToMultiPionsChannel(2,particle1, particle2);
318  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
319  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX > rChannel) {
320  isElastic = false;
321 // NN -> Eta3PiNN channel is chosen
322  INCL_DEBUG("NN interaction: Eta + three Pions channel chosen" << '\n');
323  weight = counterweight;
324  return new NNEtaToMultiPionsChannel(3,particle1, particle2);
325  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
326  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX > rChannel) {
327  isElastic = false;
328 // NN -> Eta4PiNN channel is chosen
329  INCL_DEBUG("NN interaction: Eta + four Pions channel chosen" << '\n');
330  weight = counterweight;
331  return new NNEtaToMultiPionsChannel(4,particle1, particle2);
332  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
333  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
334  + omegaProductionCX > rChannel) {
335  isElastic = false;
336 // NN -> NNOmega channel is chosen
337  INCL_DEBUG("NN interaction: Omega channel chosen" << '\n');
338  weight = counterweight;
339  return new NNToNNOmegaChannel(particle1, particle2);
340  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
341  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
342  + omegaProductionCX + omegadeltaProductionCX > rChannel) {
343  isElastic = false;
344 // NN -> N Delta Omega channel is chosen
345  INCL_DEBUG("NN interaction: Delta Omega channel chosen" << '\n');
346  weight = counterweight;
347  return new NDeltaOmegaProductionChannel(particle1, particle2);
348  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
349  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
350  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX > rChannel) {
351  isElastic = false;
352 // NN -> OmegaPiNN channel is chosen
353  INCL_DEBUG("NN interaction: Omega + one Pion channel chosen" << '\n');
354  weight = counterweight;
355  return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
356  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
357  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
358  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX > rChannel) {
359  isElastic = false;
360 // NN -> Omega2PiNN channel is chosen
361  INCL_DEBUG("NN interaction: Omega + two Pions channel chosen" << '\n');
362  weight = counterweight;
363  return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
364  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
365  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
366  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX > rChannel) {
367  isElastic = false;
368 // NN -> Omega3PiNN channel is chosen
369  INCL_DEBUG("NN interaction: Omega + three Pions channel chosen" << '\n');
370  weight = counterweight;
371  return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
372  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
373  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
374  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX > rChannel) {
375  isElastic = false;
376 // NN -> Omega4PiNN channel is chosen
377  INCL_DEBUG("NN interaction: Omega + four Pions channel chosen" << '\n');
378  weight = counterweight;
379  return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
380  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
381  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
382  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
383  + NLKProductionCX > rChannel) {
384  isElastic = false;
385  isStrangeProduction = true;
386 // NN -> NLK channel is chosen
387  INCL_DEBUG("NN interaction: NLK channel chosen" << '\n');
388  weight = bias_apply;
389  return new NNToNLKChannel(particle1, particle2);
390  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
391  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
392  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
393  + NLKProductionCX + NLKpiProductionCX > rChannel) {
394  isElastic = false;
395  isStrangeProduction = true;
396 // NN -> NLKpi channel is chosen
397  INCL_DEBUG("NN interaction: NLKpi channel chosen" << '\n');
398  weight = bias_apply;
399  return new NNToNLKpiChannel(particle1, particle2);
400  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
401  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
402  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
403  + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX > rChannel) {
404  isElastic = false;
405  isStrangeProduction = true;
406 // NN -> NLK2pi channel is chosen
407  INCL_DEBUG("NN interaction: NLK2pi channel chosen" << '\n');
408  weight = bias_apply;
409  return new NNToNLK2piChannel(particle1, particle2);
410  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
411  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
412  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
413  + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX > rChannel) {
414  isElastic = false;
415  isStrangeProduction = true;
416 // NN -> NSK channel is chosen
417  INCL_DEBUG("NN interaction: NSK channel chosen" << '\n');
418  weight = bias_apply;
419  return new NNToNSKChannel(particle1, particle2);
420  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
421  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
422  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
423  + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX > rChannel) {
424  isElastic = false;
425  isStrangeProduction = true;
426 // NN -> NSKpi channel is chosen
427  INCL_DEBUG("NN interaction: NSKpi channel chosen" << '\n');
428  weight = bias_apply;
429  return new NNToNSKpiChannel(particle1, particle2);
430  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
431  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
432  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
433  + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX > rChannel) {
434  isElastic = false;
435  isStrangeProduction = true;
436 // NN -> NSK2pi channel is chosen
437  INCL_DEBUG("NN interaction: NSK2pi channel chosen" << '\n');
438  weight = bias_apply;
439  return new NNToNSK2piChannel(particle1, particle2);
440  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
441  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
442  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
443  + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX > rChannel) {
444  isElastic = false;
445  isStrangeProduction = true;
446 // NN -> NNKKb channel is chosen
447  INCL_DEBUG("NN interaction: NNKKb channel chosen" << '\n');
448  weight = bias_apply;
449  return new NNToNNKKbChannel(particle1, particle2);
450  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
451  + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
452  + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
453  + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX> rChannel) {
454  isElastic = false;
455  isStrangeProduction = true;
456 // NN -> Missing Strangeness channel is chosen
457  INCL_DEBUG("NN interaction: Missing Strangeness channel chosen" << '\n');
458  weight = bias_apply;
459  return new NNToMissingStrangenessChannel(particle1, particle2);
460  } else {
461  INCL_WARN("inconsistency within the NN Cross Sections (sum!=inelastic)" << '\n');
462  if(NNMissingCX>0.) {
463  INCL_WARN("Returning an Missing Strangeness channel" << '\n');
464  weight = bias_apply;
465  isElastic = false;
466  isStrangeProduction = true;
467  return new NNToNNKKbChannel(particle1, particle2);
468  } else if(NNKKbProductionCX>0.) {
469  INCL_WARN("Returning an NNKKb channel" << '\n');
470  weight = bias_apply;
471  isElastic = false;
472  isStrangeProduction = true;
473  return new NNToNNKKbChannel(particle1, particle2);
474  } else if(NSK2piProductionCX>0.) {
475  INCL_WARN("Returning an NSK2pi channel" << '\n');
476  weight = bias_apply;
477  isElastic = false;
478  isStrangeProduction = true;
479  return new NNToNSK2piChannel(particle1, particle2);
480  } else if(NSKpiProductionCX>0.) {
481  INCL_WARN("Returning an NSKpi channel" << '\n');
482  weight = bias_apply;
483  isElastic = false;
484  isStrangeProduction = true;
485  return new NNToNSKpiChannel(particle1, particle2);
486  } else if(NSKProductionCX>0.) {
487  INCL_WARN("Returning an NSK channel" << '\n');
488  weight = bias_apply;
489  isElastic = false;
490  isStrangeProduction = true;
491  return new NNToNSKChannel(particle1, particle2);
492  } else if(NLK2piProductionCX>0.) {
493  INCL_WARN("Returning an NLK2pi channel" << '\n');
494  weight = bias_apply;
495  isElastic = false;
496  isStrangeProduction = true;
497  return new NNToNLK2piChannel(particle1, particle2);
498  } else if(NLKpiProductionCX>0.) {
499  INCL_WARN("Returning an NLKpi channel" << '\n');
500  weight = bias_apply;
501  isElastic = false;
502  isStrangeProduction = true;
503  return new NNToNLKpiChannel(particle1, particle2);
504  } else if(NLKProductionCX>0.) {
505  INCL_WARN("Returning an NLK channel" << '\n');
506  weight = bias_apply;
507  isElastic = false;
508  isStrangeProduction = true;
509  return new NNToNLKChannel(particle1, particle2);
510  } else if(omegafourPiProductionCX>0.) {
511  INCL_WARN("Returning an Omega + four Pions channel" << '\n');
512  weight = counterweight;
513  isElastic = false;
514  return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
515  } else if(omegathreePiProductionCX>0.) {
516  INCL_WARN("Returning an Omega + three Pions channel" << '\n');
517  weight = counterweight;
518  isElastic = false;
519  return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
520  } else if(omegatwoPiProductionCX>0.) {
521  INCL_WARN("Returning an Omega + two Pions channel" << '\n');
522  weight = counterweight;
523  isElastic = false;
524  return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
525  } else if(omegaonePiProductionCX>0.) {
526  INCL_WARN("Returning an Omega + one Pion channel" << '\n');
527  weight = counterweight;
528  isElastic = false;
529  return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
530  } else if(omegadeltaProductionCX>0.) {
531  INCL_WARN("Returning an Omega + Delta channel" << '\n');
532  weight = counterweight;
533  isElastic = false;
534  return new NDeltaOmegaProductionChannel(particle1, particle2);
535  } else if(omegaProductionCX>0.) {
536  INCL_WARN("Returning an Omega channel" << '\n');
537  weight = counterweight;
538  isElastic = false;
539  return new NNToNNOmegaChannel(particle1, particle2);
540  } else if(etafourPiProductionCX>0.) {
541  INCL_WARN("Returning an Eta + four Pions channel" << '\n');
542  weight = counterweight;
543  isElastic = false;
544  return new NNEtaToMultiPionsChannel(4,particle1, particle2);
545  } else if(etathreePiProductionCX>0.) {
546  INCL_WARN("Returning an Eta + threev channel" << '\n');
547  weight = counterweight;
548  isElastic = false;
549  return new NNEtaToMultiPionsChannel(3,particle1, particle2);
550  } else if(etatwoPiProductionCX>0.) {
551  INCL_WARN("Returning an Eta + two Pions channel" << '\n');
552  weight = counterweight;
553  isElastic = false;
554  return new NNEtaToMultiPionsChannel(2,particle1, particle2);
555  } else if(etaonePiProductionCX>0.) {
556  INCL_WARN("Returning an Eta + one Pion channel" << '\n');
557  weight = counterweight;
558  isElastic = false;
559  return new NNEtaToMultiPionsChannel(1,particle1, particle2);
560  } else if(etadeltaProductionCX>0.) {
561  INCL_WARN("Returning an Eta + Delta channel" << '\n');
562  weight = counterweight;
563  isElastic = false;
564  return new NDeltaEtaProductionChannel(particle1, particle2);
565  } else if(etaProductionCX>0.) {
566  INCL_WARN("Returning an Eta channel" << '\n');
567  weight = counterweight;
568  isElastic = false;
569  return new NNToNNEtaChannel(particle1, particle2);
570  } else if(fourPiProductionCX>0.) {
571  INCL_WARN("Returning a 4pi channel" << '\n');
572  weight = counterweight;
573  isElastic = false;
574  return new NNToMultiPionsChannel(4,particle1, particle2);
575  } else if(threePiProductionCX>0.) {
576  INCL_WARN("Returning a 3pi channel" << '\n');
577  weight = counterweight;
578  isElastic = false;
579  return new NNToMultiPionsChannel(3,particle1, particle2);
580  } else if(twoPiProductionCX>0.) {
581  INCL_WARN("Returning a 2pi channel" << '\n');
582  weight = counterweight;
583  isElastic = false;
584  return new NNToMultiPionsChannel(2,particle1, particle2);
585  } else if(onePiProductionCX>0.) {
586  INCL_WARN("Returning a 1pi channel" << '\n');
587  weight = counterweight;
588  isElastic = false;
589  return new NNToMultiPionsChannel(1,particle1, particle2);
590  } else if(deltaProductionCX>0.) {
591  INCL_WARN("Returning a delta-production channel" << '\n');
592  weight = counterweight;
593  isElastic = false;
594  return new DeltaProductionChannel(particle1, particle2);
595  } else {
596  INCL_WARN("Returning an elastic channel" << '\n');
597  weight = counterweight;
598  isElastic = true;
599  return new ElasticChannel(particle1, particle2);
600  }
601  }
602 
604  }
605  else if((particle1->isNucleon() && particle2->isDelta()) ||
606  (particle1->isDelta() && particle2->isNucleon())) {
607 
608  G4double NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
609  G4double NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
610  G4double DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
611  G4double DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
612  G4double NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
613 
614  const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NDeltaToNN(particle1, particle2);
615  const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX)/bias_apply;
616 
617  G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
618 
619  if(counterweight < 0.5){
620  counterweight = 0.5;
621  bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
622 
623  NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
624  NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
625  DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
626  DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
627  NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
628  }
629 
630  G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
631  G4double recombinationCX = CrossSections::NDeltaToNN(particle1, particle2)*counterweight;
632 
633  const G4double rChannel=Random::shoot() * (StrangenessProdCX + UnStrangeProdCX);
634 
635  if(elasticCX > rChannel) {
636  isElastic = true;
637 // Elastic N Delta channel
638  INCL_DEBUG("NDelta interaction: elastic channel chosen" << '\n');
639  weight = counterweight;
640  return new ElasticChannel(particle1, particle2);
641  } else if (elasticCX + recombinationCX > rChannel){
642  isElastic = false;
643 // Recombination
644 // NDelta -> NN channel is chosen
645  INCL_DEBUG("NDelta interaction: recombination channel chosen" << '\n');
646  weight = counterweight;
647  return new RecombinationChannel(particle1, particle2);
648  } else if (elasticCX + recombinationCX + NLKProductionCX > rChannel){
649  isElastic = false;
650  isStrangeProduction = true;
651 // NDelta -> NLK channel is chosen
652  INCL_DEBUG("NDelta interaction: NLK channel chosen" << '\n');
653  weight = bias_apply;
654  return new NDeltaToNLKChannel(particle1, particle2);
655  } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX > rChannel){
656  isElastic = false;
657  isStrangeProduction = true;
658 // NDelta -> NSK channel is chosen
659  INCL_DEBUG("NDelta interaction: NSK channel chosen" << '\n');
660  weight = bias_apply;
661  return new NDeltaToNSKChannel(particle1, particle2);
662  } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX > rChannel){
663  isElastic = false;
664  isStrangeProduction = true;
665 // NDelta -> DeltaLK channel is chosen
666  INCL_DEBUG("NDelta interaction: DeltaLK channel chosen" << '\n');
667  weight = bias_apply;
668  return new NDeltaToDeltaLKChannel(particle1, particle2);
669  } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX > rChannel){
670  isElastic = false;
671  isStrangeProduction = true;
672 // NDelta -> DeltaSK channel is chosen
673  INCL_DEBUG("NDelta interaction: DeltaSK channel chosen" << '\n');
674  weight = bias_apply;
675  return new NDeltaToDeltaSKChannel(particle1, particle2);
676  } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX > rChannel){
677  isElastic = false;
678  isStrangeProduction = true;
679 // NDelta -> NNKKb channel is chosen
680  INCL_DEBUG("NDelta interaction: NNKKb channel chosen" << '\n');
681  weight = bias_apply;
682  return new NDeltaToNNKKbChannel(particle1, particle2);
683  }
684  else{
685  INCL_ERROR("rChannel > (StrangenessProdCX + UnStrangeProdCX) in NDelta interaction: return an elastic channel" << '\n');
686  weight = counterweight;
687  isElastic = true;
688  return new ElasticChannel(particle1, particle2);
689  }
690 
692  } else if(particle1->isDelta() && particle2->isDelta()) {
693  isElastic = true;
694  INCL_DEBUG("DeltaDelta interaction: elastic channel chosen" << '\n');
695  return new ElasticChannel(particle1, particle2);
696 
698  } else if(isPiN) {
699 
700  G4double LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
701  G4double SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
702  G4double LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
703  G4double SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
704  G4double LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
705  G4double SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
706  G4double NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
707  G4double MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
708 
709  const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::piNToDelta(particle1, particle2)
710  + CrossSections::piNToxPiN(2,particle1, particle2) + CrossSections::piNToxPiN(3,particle1, particle2) + CrossSections::piNToxPiN(4,particle1, particle2)
711  + CrossSections::piNToEtaN(particle1, particle2) + CrossSections::piNToOmegaN(particle1, particle2);
712  const G4double StrangenessProdCX = (LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX)/bias_apply;
713 
714  G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
715 
716  if(counterweight < 0.5) {
717  counterweight = 0.5;
718  bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
719  LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
720  SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
721  LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
722  SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
723  LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
724  SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
725  NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
726  MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
727  }
728 
729 
730  const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
731  const G4double deltaProductionCX = CrossSections::piNToDelta(particle1, particle2)*counterweight;
732  const G4double onePiProductionCX = CrossSections::piNToxPiN(2,particle1, particle2)*counterweight;
733  const G4double twoPiProductionCX = CrossSections::piNToxPiN(3,particle1, particle2)*counterweight;
734  const G4double threePiProductionCX = CrossSections::piNToxPiN(4,particle1, particle2)*counterweight;
735  const G4double etaProductionCX = CrossSections::piNToEtaN(particle1, particle2)*counterweight;
736  const G4double omegaProductionCX = CrossSections::piNToOmegaN(particle1, particle2)*counterweight;
737 
738  const G4double totCX=CrossSections::total(particle1, particle2);
739 
740 // assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-etaProductionCX-omegaProductionCX-LKProdCX-SKProdCX-LKpiProdCX-SKpiProdCX-LK2piProdCX-SK2piProdCX-NKKbProdCX-MissingCX) < 0.15);
741 
742  const G4double rChannel=Random::shoot() * totCX;
743 
744  if(elasticCX > rChannel) {
745  isElastic = true;
746 // Elastic PiN channel
747  INCL_DEBUG("PiN interaction: elastic channel chosen" << '\n');
748  weight = counterweight;
749  return new PiNElasticChannel(particle1, particle2);
750  } else if(elasticCX + deltaProductionCX > rChannel) {
751  isElastic = false;
752 // PiN -> Delta channel is chosen
753  INCL_DEBUG("PiN interaction: Delta channel chosen" << '\n');
754  weight = counterweight;
755  return new PiNToDeltaChannel(particle1, particle2);
756  } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
757  isElastic = false;
758 // PiN -> PiNPi channel is chosen
759  INCL_DEBUG("PiN interaction: one Pion channel chosen" << '\n');
760  weight = counterweight;
761  return new PiNToMultiPionsChannel(2,particle1, particle2);
762  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
763  isElastic = false;
764 // PiN -> PiN2Pi channel is chosen
765  INCL_DEBUG("PiN interaction: two Pions channel chosen" << '\n');
766  weight = counterweight;
767  return new PiNToMultiPionsChannel(3,particle1, particle2);
768  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
769  isElastic = false;
770 // PiN -> PiN3Pi channel is chosen
771  INCL_DEBUG("PiN interaction: three Pions channel chosen" << '\n');
772  weight = counterweight;
773  return new PiNToMultiPionsChannel(4,particle1, particle2);
774  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX > rChannel) {
775  isElastic = false;
776 // PiN -> EtaN channel is chosen
777  INCL_DEBUG("PiN interaction: Eta channel chosen" << '\n');
778  weight = counterweight;
779  return new PiNToEtaChannel(particle1, particle2);
780  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX > rChannel) {
781  isElastic = false;
782 // PiN -> OmegaN channel is chosen
783  INCL_DEBUG("PiN interaction: Omega channel chosen" << '\n');
784  weight = counterweight;
785  return new PiNToOmegaChannel(particle1, particle2);
786  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
787  + LKProdCX > rChannel) {
788  isElastic = false;
789  isStrangeProduction = true;
790 // PiN -> LK channel is chosen
791  INCL_DEBUG("PiN interaction: LK channel chosen" << '\n');
792  weight = bias_apply;
793  return new NpiToLKChannel(particle1, particle2);
794  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
795  + LKProdCX + SKProdCX > rChannel) {
796  isElastic = false;
797  isStrangeProduction = true;
798 // PiN -> SK channel is chosen
799  INCL_DEBUG("PiN interaction: SK channel chosen" << '\n');
800  weight = bias_apply;
801  return new NpiToSKChannel(particle1, particle2);
802  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
803  + LKProdCX + SKProdCX + LKpiProdCX > rChannel) {
804  isElastic = false;
805  isStrangeProduction = true;
806 // PiN -> LKpi channel is chosen
807  INCL_DEBUG("PiN interaction: LKpi channel chosen" << '\n');
808  weight = bias_apply;
809  return new NpiToLKpiChannel(particle1, particle2);
810  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
811  + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX > rChannel) {
812  isElastic = false;
813  isStrangeProduction = true;
814 // PiN -> SKpi channel is chosen
815  INCL_DEBUG("PiN interaction: SKpi channel chosen" << '\n');
816  weight = bias_apply;
817  return new NpiToSKpiChannel(particle1, particle2);
818  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
819  + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX > rChannel) {
820  isElastic = false;
821  isStrangeProduction = true;
822 // PiN -> LK2pi channel is chosen
823  INCL_DEBUG("PiN interaction: LK2pi channel chosen" << '\n');
824  weight = bias_apply;
825  return new NpiToLK2piChannel(particle1, particle2);
826  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
827  + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX > rChannel) {
828  isElastic = false;
829  isStrangeProduction = true;
830 // PiN -> SK2pi channel is chosen
831  INCL_DEBUG("PiN interaction: SK2pi channel chosen" << '\n');
832  weight = bias_apply;
833  return new NpiToSK2piChannel(particle1, particle2);
834  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
835  + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX > rChannel) {
836  isElastic = false;
837  isStrangeProduction = true;
838 // PiN -> NKKb channel is chosen
839  INCL_DEBUG("PiN interaction: NKKb channel chosen" << '\n');
840  weight = bias_apply;
841  return new NpiToNKKbChannel(particle1, particle2);
842  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
843  + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX> rChannel) {
844  isElastic = false;
845  isStrangeProduction = true;
846 // PiN -> Missinge Strangeness channel is chosen
847  INCL_DEBUG("PiN interaction: Missinge Strangeness channel chosen" << '\n');
848  weight = bias_apply;
849  return new NpiToMissingStrangenessChannel(particle1, particle2);
850  }
851  else {
852  INCL_WARN("inconsistency within the PiN Cross Sections (sum!=inelastic)" << '\n');
853  if(MissingCX>0.) {
854  INCL_WARN("Returning a Missinge Strangeness channel" << '\n');
855  weight = bias_apply;
856  isElastic = false;
857  isStrangeProduction = true;
858  return new NpiToMissingStrangenessChannel(particle1, particle2);
859  } else if(NKKbProdCX>0.) {
860  INCL_WARN("Returning a NKKb channel" << '\n');
861  weight = bias_apply;
862  isElastic = false;
863  isStrangeProduction = true;
864  return new NpiToNKKbChannel(particle1, particle2);
865  } else if(SK2piProdCX>0.) {
866  INCL_WARN("Returning a SK2pi channel" << '\n');
867  weight = bias_apply;
868  isElastic = false;
869  isStrangeProduction = true;
870  return new NpiToSK2piChannel(particle1, particle2);
871  } else if(LK2piProdCX>0.) {
872  INCL_WARN("Returning a LK2pi channel" << '\n');
873  weight = bias_apply;
874  isElastic = false;
875  isStrangeProduction = true;
876  return new NpiToLK2piChannel(particle1, particle2);
877  } else if(SKpiProdCX>0.) {
878  INCL_WARN("Returning a SKpi channel" << '\n');
879  weight = bias_apply;
880  isElastic = false;
881  isStrangeProduction = true;
882  return new NpiToSKpiChannel(particle1, particle2);
883  } else if(LKpiProdCX>0.) {
884  INCL_WARN("Returning a LKpi channel" << '\n');
885  weight = bias_apply;
886  isElastic = false;
887  isStrangeProduction = true;
888  return new NpiToLKpiChannel(particle1, particle2);
889  } else if(SKProdCX>0.) {
890  INCL_WARN("Returning a SK channel" << '\n');
891  weight = bias_apply;
892  isElastic = false;
893  isStrangeProduction = true;
894  return new NpiToSKChannel(particle1, particle2);
895  } else if(LKProdCX>0.) {
896  INCL_WARN("Returning a LK channel" << '\n');
897  weight = bias_apply;
898  isElastic = false;
899  isStrangeProduction = true;
900  return new NpiToLKChannel(particle1, particle2);
901  } else if(omegaProductionCX>0.) {
902  INCL_WARN("Returning a Omega channel" << '\n');
903  weight = counterweight;
904  isElastic = false;
905  return new PiNToOmegaChannel(particle1, particle2);
906  } else if(etaProductionCX>0.) {
907  INCL_WARN("Returning a Eta channel" << '\n');
908  weight = counterweight;
909  isElastic = false;
910  return new PiNToEtaChannel(particle1, particle2);
911  } else if(threePiProductionCX>0.) {
912  INCL_WARN("Returning a 3pi channel" << '\n');
913  weight = counterweight;
914  isElastic = false;
915  return new PiNToMultiPionsChannel(4,particle1, particle2);
916  } else if(twoPiProductionCX>0.) {
917  INCL_WARN("Returning a 2pi channel" << '\n');
918  weight = counterweight;
919  isElastic = false;
920  return new PiNToMultiPionsChannel(3,particle1, particle2);
921  } else if(onePiProductionCX>0.) {
922  INCL_WARN("Returning a 1pi channel" << '\n');
923  weight = counterweight;
924  isElastic = false;
925  return new PiNToMultiPionsChannel(2,particle1, particle2);
926  } else if(deltaProductionCX>0.) {
927  INCL_WARN("Returning a delta-production channel" << '\n');
928  weight = counterweight;
929  isElastic = false;
930  return new PiNToDeltaChannel(particle1, particle2);
931  } else {
932  INCL_WARN("Returning an elastic channel" << '\n');
933  weight = counterweight;
934  isElastic = true;
935  return new PiNElasticChannel(particle1, particle2);
936  }
937  }
938  } else if ((particle1->isNucleon() && particle2->isEta()) || (particle2->isNucleon() && particle1->isEta())) {
940 
941  const G4double elasticCX = CrossSections::elastic(particle1, particle2);
942  const G4double onePiProductionCX = CrossSections::etaNToPiN(particle1, particle2);
943  const G4double twoPiProductionCX = CrossSections::etaNToPiPiN(particle1, particle2);
944  const G4double totCX=CrossSections::total(particle1, particle2);
945 // assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
946 
947  const G4double rChannel=Random::shoot() * totCX;
948 
949  if(elasticCX > rChannel) {
950 // Elastic EtaN channel
951  isElastic = true;
952  INCL_DEBUG("EtaN interaction: elastic channel chosen" << '\n');
953  return new EtaNElasticChannel(particle1, particle2);
954  } else if(elasticCX + onePiProductionCX > rChannel) {
955  isElastic = false;
956 // EtaN -> EtaPiN channel is chosen
957  INCL_DEBUG("EtaN interaction: PiN channel chosen" << '\n');
958  return new EtaNToPiNChannel(particle1, particle2);
959  } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
960  isElastic = false;
961 // EtaN -> EtaPiPiN channel is chosen
962  INCL_DEBUG("EtaN interaction: PiPiN channel chosen" << '\n');
963  return new EtaNToPiPiNChannel(particle1, particle2);
964  }
965 
966  else {
967  INCL_WARN("inconsistency within the EtaN Cross Sections (sum!=inelastic)" << '\n');
968  if(twoPiProductionCX>0.) {
969  INCL_WARN("Returning a PiPiN channel" << '\n');
970  isElastic = false;
971  return new EtaNToPiPiNChannel(particle1, particle2);
972  } else if(onePiProductionCX>0.) {
973  INCL_WARN("Returning a PiN channel" << '\n');
974  isElastic = false;
975  return new EtaNToPiNChannel(particle1, particle2);
976  } else {
977  INCL_WARN("Returning an elastic channel" << '\n');
978  isElastic = true;
979  return new EtaNElasticChannel(particle1, particle2);
980  }
981  }
982 
983  } else if ((particle1->isNucleon() && particle2->isOmega()) || (particle2->isNucleon() && particle1->isOmega())) {
985 
986  const G4double elasticCX = CrossSections::elastic(particle1, particle2);
987  const G4double onePiProductionCX = CrossSections::omegaNToPiN(particle1, particle2);
988  const G4double twoPiProductionCX = CrossSections::omegaNToPiPiN(particle1, particle2);
989  const G4double totCX=CrossSections::total(particle1, particle2);
990 // assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
991 
992  const G4double rChannel=Random::shoot() * totCX;
993 
994  if(elasticCX > rChannel) {
995 // Elastic OmegaN channel
996  isElastic = true;
997  INCL_DEBUG("OmegaN interaction: elastic channel chosen" << '\n');
998  return new OmegaNElasticChannel(particle1, particle2);
999  } else if(elasticCX + onePiProductionCX > rChannel) {
1000  isElastic = false;
1001 // OmegaN -> PiN channel is chosen
1002  INCL_DEBUG("OmegaN interaction: PiN channel chosen" << '\n');
1003  return new OmegaNToPiNChannel(particle1, particle2);
1004  } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
1005  isElastic = false;
1006 // OmegaN -> PiPiN channel is chosen
1007  INCL_DEBUG("OmegaN interaction: PiPiN channel chosen" << '\n');
1008  return new OmegaNToPiPiNChannel(particle1, particle2);
1009  }
1010  else {
1011  INCL_WARN("inconsistency within the OmegaN Cross Sections (sum!=inelastic)" << '\n');
1012  if(twoPiProductionCX>0.) {
1013  INCL_WARN("Returning a PiPiN channel" << '\n');
1014  isElastic = false;
1015  return new OmegaNToPiPiNChannel(particle1, particle2);
1016  } else if(onePiProductionCX>0.) {
1017  INCL_WARN("Returning a PiN channel" << '\n');
1018  isElastic = false;
1019  return new OmegaNToPiNChannel(particle1, particle2);
1020  } else {
1021  INCL_WARN("Returning an elastic channel" << '\n');
1022  isElastic = true;
1023  return new OmegaNElasticChannel(particle1, particle2);
1024  }
1025  }
1026  } else if ((particle1->isNucleon() && particle2->isKaon()) || (particle2->isNucleon() && particle1->isKaon())) {
1028  const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1029  const G4double quasielasticCX = CrossSections::NKToNK(particle1,particle2);
1030  const G4double NKToNKpiCX = CrossSections::NKToNKpi(particle1,particle2);
1031  const G4double NKToNK2piCX = CrossSections::NKToNK2pi(particle1,particle2);
1032  const G4double totCX=CrossSections::total(particle1, particle2);
1033 // assert(std::fabs(totCX-elasticCX-quasielasticCX-NKToNKpiCX-NKToNK2piCX)<0.1);
1034 
1035  const G4double rChannel=Random::shoot() * totCX;
1036  if(elasticCX > rChannel){
1037 // Elastic KN channel is chosen
1038  isElastic = true;
1039  INCL_DEBUG("KN interaction: elastic channel chosen" << '\n');
1040  return new NKElasticChannel(particle1, particle2);
1041  } else if(elasticCX + quasielasticCX > rChannel){
1042 // Quasi-elastic KN channel is chosen
1043  isElastic = false; // true ??
1044  INCL_DEBUG("KN interaction: quasi-elastic channel chosen" << '\n');
1045  return new NKToNKChannel(particle1, particle2);
1046  } else if(elasticCX + quasielasticCX + NKToNKpiCX > rChannel){
1047 // KN -> NKpi channel is chosen
1048  isElastic = false;
1049  INCL_DEBUG("KN interaction: NKpi channel chosen" << '\n');
1050  return new NKToNKpiChannel(particle1, particle2);
1051  } else if(elasticCX + quasielasticCX + NKToNKpiCX + NKToNK2piCX > rChannel){
1052 // KN -> NK2pi channel is chosen
1053  isElastic = false;
1054  INCL_DEBUG("KN interaction: NK2pi channel chosen" << '\n');
1055  return new NKToNK2piChannel(particle1, particle2);
1056  } else {
1057  INCL_WARN("inconsistency within the KN Cross Sections (sum!=inelastic)" << '\n');
1058  if(NKToNK2piCX>0.) {
1059  INCL_WARN("Returning a NKToNK2pi channel" << '\n');
1060  isElastic = false;
1061  return new NKToNK2piChannel(particle1, particle2);
1062  } else if(NKToNKpiCX>0.) {
1063  INCL_WARN("Returning a NKToNKpi channel" << '\n');
1064  isElastic = false;
1065  return new NKToNKpiChannel(particle1, particle2);
1066  } else if(quasielasticCX>0.) {
1067  INCL_WARN("Returning a quasi-elastic channel" << '\n');
1068  isElastic = false; // true ??
1069  return new NKToNKChannel(particle1, particle2);
1070  } else {
1071  INCL_WARN("Returning an elastic channel" << '\n');
1072  isElastic = true;
1073  return new NKElasticChannel(particle1, particle2);
1074  }
1075  }
1076  } else if ((particle1->isNucleon() && particle2->isAntiKaon()) || (particle2->isNucleon() && particle1->isAntiKaon())) {
1078  const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1079  const G4double quasielasticCX = CrossSections::NKbToNKb(particle1,particle2);
1080  const G4double NKbToNKbpiCX = CrossSections::NKbToNKbpi(particle1,particle2);
1081  const G4double NKbToNKb2piCX = CrossSections::NKbToNKb2pi(particle1,particle2);
1082  const G4double NKbToLpiCX = CrossSections::NKbToLpi(particle1,particle2);
1083  const G4double NKbToL2piCX = CrossSections::NKbToL2pi(particle1,particle2);
1084  const G4double NKbToSpiCX = CrossSections::NKbToSpi(particle1,particle2);
1085  const G4double NKbToS2piCX = CrossSections::NKbToS2pi(particle1,particle2);
1086  const G4double totCX=CrossSections::total(particle1, particle2);
1087 // assert(std::fabs(totCX-elasticCX-quasielasticCX-NKbToNKbpiCX-NKbToNKb2piCX-NKbToLpiCX-NKbToL2piCX-NKbToSpiCX-NKbToS2piCX)<0.1);
1088 
1089  const G4double rChannel=Random::shoot() * totCX;
1090  if(elasticCX > rChannel){
1091 // Elastic KbN channel is chosen
1092  isElastic = true;
1093  INCL_DEBUG("KbN interaction: elastic channel chosen" << '\n');
1094  return new NKbElasticChannel(particle1, particle2);
1095  } else if(elasticCX + quasielasticCX > rChannel){
1096 // Quasi-elastic KbN channel is chosen
1097  isElastic = false; // true ??
1098  INCL_DEBUG("KbN interaction: quasi-elastic channel chosen" << '\n');
1099  return new NKbToNKbChannel(particle1, particle2);
1100  } else if(elasticCX + quasielasticCX + NKbToNKbpiCX > rChannel){
1101 // KbN -> NKbpi channel is chosen
1102  isElastic = false;
1103  INCL_DEBUG("KbN interaction: NKbpi channel chosen" << '\n');
1104  return new NKbToNKbpiChannel(particle1, particle2);
1105  } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX > rChannel){
1106 // KbN -> NKb2pi channel is chosen
1107  isElastic = false;
1108  INCL_DEBUG("KbN interaction: NKb2pi channel chosen" << '\n');
1109  return new NKbToNKb2piChannel(particle1, particle2);
1110  } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX > rChannel){
1111 // KbN -> Lpi channel is chosen
1112  isElastic = false;
1113  INCL_DEBUG("KbN interaction: Lpi channel chosen" << '\n');
1114  return new NKbToLpiChannel(particle1, particle2);
1115  } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX > rChannel){
1116 // KbN -> L2pi channel is chosen
1117  isElastic = false;
1118  INCL_DEBUG("KbN interaction: L2pi channel chosen" << '\n');
1119  return new NKbToL2piChannel(particle1, particle2);
1120  } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX > rChannel){
1121 // KbN -> Spi channel is chosen
1122  isElastic = false;
1123  INCL_DEBUG("KbN interaction: Spi channel chosen" << '\n');
1124  return new NKbToSpiChannel(particle1, particle2);
1125  } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX + NKbToS2piCX > rChannel){
1126 // KbN -> S2pi channel is chosen
1127  isElastic = false;
1128  INCL_DEBUG("KbN interaction: S2pi channel chosen" << '\n');
1129  return new NKbToS2piChannel(particle1, particle2);
1130  } else {
1131  INCL_WARN("inconsistency within the KbN Cross Sections (sum!=inelastic)" << '\n');
1132  if(NKbToS2piCX>0.) {
1133  INCL_WARN("Returning a NKbToS2pi channel" << '\n');
1134  isElastic = false;
1135  return new NKbToS2piChannel(particle1, particle2);
1136  } else if(NKbToSpiCX>0.) {
1137  INCL_WARN("Returning a NKbToSpi channel" << '\n');
1138  isElastic = false;
1139  return new NKbToSpiChannel(particle1, particle2);
1140  } else if(NKbToL2piCX>0.) {
1141  INCL_WARN("Returning a NKbToL2pi channel" << '\n');
1142  isElastic = false;
1143  return new NKbToL2piChannel(particle1, particle2);
1144  } else if(NKbToLpiCX>0.) {
1145  INCL_WARN("Returning a NKbToLpi channel" << '\n');
1146  isElastic = false;
1147  return new NKbToLpiChannel(particle1, particle2);
1148  } else if(NKbToNKb2piCX>0.) {
1149  INCL_WARN("Returning a NKbToNKb2pi channel" << '\n');
1150  isElastic = false;
1151  return new NKbToNKb2piChannel(particle1, particle2);
1152  } else if(NKbToNKbpiCX>0.) {
1153  INCL_WARN("Returning a NKbToNKbpi channel" << '\n');
1154  isElastic = false;
1155  return new NKbToNKbpiChannel(particle1, particle2);
1156  } else if(quasielasticCX>0.) {
1157  INCL_WARN("Returning a quasi-elastic channel" << '\n');
1158  isElastic = false; // true ??
1159  return new NKbToNKbChannel(particle1, particle2);
1160  } else {
1161  INCL_WARN("Returning an elastic channel" << '\n');
1162  isElastic = true;
1163  return new NKbElasticChannel(particle1, particle2);
1164  }
1165  }
1166  } else if ((particle1->isNucleon() && particle2->isLambda()) || (particle2->isNucleon() && particle1->isLambda())) {
1168  const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1169  const G4double NLToNSCX = CrossSections::NLToNS(particle1,particle2);
1170  const G4double totCX=CrossSections::total(particle1, particle2);
1171 // assert(std::fabs(totCX-elasticCX-NLToNSCX)<0.1);
1172 
1173  const G4double rChannel=Random::shoot() * totCX;
1174  if(elasticCX > rChannel){
1175 // Elastic NLambda channel is chosen
1176  isElastic = true;
1177  INCL_DEBUG("NLambda interaction: elastic channel chosen" << '\n');
1178  return new NYElasticChannel(particle1, particle2);
1179  } else if(elasticCX + NLToNSCX > rChannel){
1180 // Quasi-elastic NLambda channel is chosen
1181  isElastic = false; // true ??
1182  INCL_DEBUG("NLambda interaction: quasi-elastic channel chosen" << '\n');
1183  return new NLToNSChannel(particle1, particle2);
1184  } else {
1185  INCL_WARN("inconsistency within the NLambda Cross Sections (sum!=inelastic)" << '\n');
1186  if(NLToNSCX>0.) {
1187  INCL_WARN("Returning a quasi-elastic channel" << '\n');
1188  isElastic = false; // true ??
1189  return new NLToNSChannel(particle1, particle2);
1190  } else {
1191  INCL_WARN("Returning an elastic channel" << '\n');
1192  isElastic = true;
1193  return new NYElasticChannel(particle1, particle2);
1194  }
1195  }
1196  } else if ((particle1->isNucleon() && particle2->isSigma()) || (particle2->isNucleon() && particle1->isSigma())) {
1198  const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1199  const G4double NSToNLCX = CrossSections::NSToNL(particle1,particle2);
1200  const G4double NSToNSCX = CrossSections::NSToNS(particle1,particle2);
1201  const G4double totCX=CrossSections::total(particle1, particle2);
1202 // assert(std::fabs(totCX-elasticCX-NSToNLCX-NSToNSCX)<0.1);
1203 
1204  const G4double rChannel=Random::shoot() * totCX;
1205  if(elasticCX > rChannel){
1206 // Elastic NSigma channel is chosen
1207  isElastic = true;
1208  INCL_DEBUG("NSigma interaction: elastic channel chosen" << '\n');
1209  return new NYElasticChannel(particle1, particle2);
1210  } else if(elasticCX + NSToNLCX > rChannel){
1211 // NSigma -> NLambda channel is chosen
1212  isElastic = false; // true ??
1213  INCL_DEBUG("NSigma interaction: NLambda channel chosen" << '\n');
1214  return new NSToNLChannel(particle1, particle2);
1215  } else if(elasticCX + NSToNLCX + NSToNSCX > rChannel){
1216 // NSigma -> NSigma quasi-elastic channel is chosen
1217  isElastic = false; // true ??
1218  INCL_DEBUG("NSigma interaction: NSigma quasi-elastic channel chosen" << '\n');
1219  return new NSToNSChannel(particle1, particle2);
1220  } else {
1221  INCL_WARN("inconsistency within the NSigma Cross Sections (sum!=inelastic)" << '\n');
1222  if(NSToNSCX>0.) {
1223  INCL_WARN("Returning a quasi-elastic channel" << '\n');
1224  isElastic = false; // true ??
1225  return new NSToNSChannel(particle1, particle2);
1226  } else if(NSToNLCX>0.) {
1227  INCL_WARN("Returning a NLambda channel" << '\n');
1228  isElastic = false; // true ??
1229  return new NSToNLChannel(particle1, particle2);
1230  } else {
1231  INCL_WARN("Returning an elastic channel" << '\n');
1232  isElastic = true;
1233  return new NYElasticChannel(particle1, particle2);
1234  }
1235  }
1236  }
1237 
1238  else {
1239  INCL_DEBUG("BinaryCollisionAvatar can only handle nucleons (for the moment)."
1240  << '\n'
1241  << particle1->print()
1242  << '\n'
1243  << particle2->print()
1244  << '\n');
1245  InteractionAvatar::restoreParticles();
1246  return NULL;
1247  }
1248  }
1249 
1250  void BinaryCollisionAvatar::preInteraction() {
1251  isParticle1Spectator = particle1->isTargetSpectator();
1252  isParticle2Spectator = particle2->isTargetSpectator();
1253  InteractionAvatar::preInteraction();
1254  }
1255 
1256  void BinaryCollisionAvatar::postInteraction(FinalState *fs) {
1257  // Call the postInteraction method of the parent class
1258  // (provides Pauli blocking and enforces energy conservation)
1259  InteractionAvatar::postInteraction(fs);
1260 
1261  switch(fs->getValidity()) {
1262  case PauliBlockedFS:
1263  theNucleus->getStore()->getBook().incrementBlockedCollisions();
1264  break;
1265  case NoEnergyConservationFS:
1266  case ParticleBelowFermiFS:
1267  case ParticleBelowZeroFS:
1268  break;
1269  case ValidFS:
1270  Book &theBook = theNucleus->getStore()->getBook();
1271  theBook.incrementAcceptedCollisions();
1272 
1273  if(theBook.getAcceptedCollisions() == 1) {
1274  // Store time and cross section of the first collision
1275  G4double t = theBook.getCurrentTime();
1276  theBook.setFirstCollisionTime(t);
1277  theBook.setFirstCollisionXSec(oldXSec);
1278  // Increase the number of Kaon by 1
1279  if(isStrangeProduction) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()+1);
1280  // Store position and momentum of the spectator on the first
1281  // collision
1282  if((isParticle1Spectator && isParticle2Spectator) || (!isParticle1Spectator && !isParticle2Spectator)) {
1283  INCL_ERROR("First collision must be within a target spectator and a non-target spectator");
1284  }
1285  if(isParticle1Spectator) {
1286  theBook.setFirstCollisionSpectatorPosition(backupParticle1->getPosition().mag());
1287  theBook.setFirstCollisionSpectatorMomentum(backupParticle1->getMomentum().mag());
1288  } else {
1289  theBook.setFirstCollisionSpectatorPosition(backupParticle2->getPosition().mag());
1290  theBook.setFirstCollisionSpectatorMomentum(backupParticle2->getMomentum().mag());
1291  }
1292 
1293  // Store the elasticity of the first collision
1294  theBook.setFirstCollisionIsElastic(isElastic);
1295  }
1296  }
1297  return;
1298  }
1299 
1300  std::string BinaryCollisionAvatar::dump() const {
1301  std::stringstream ss;
1302  ss << "(avatar " << theTime <<" 'nn-collision" << '\n'
1303  << "(list " << '\n'
1304  << particle1->dump()
1305  << particle2->dump()
1306  << "))" << '\n';
1307  return ss.str();
1308  }
1309 
1310 }