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