Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4INCLBinaryCollisionAvatar.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 /*
39  * G4INCLBinaryCollisionAvatar.cc
40  *
41  * \date Jun 5, 2009
42  * \author Pekka Kaitaniemi
43  */
44 
45 #include "G4INCLBinaryCollisionAvatar.hh"
46 #include "G4INCLElasticChannel.hh"
47 #include "G4INCLRecombinationChannel.hh"
48 #include "G4INCLDeltaProductionChannel.hh"
49 #include "G4INCLNNToMultiPionsChannel.hh"
50 #include "G4INCLNNToNNEtaChannel.hh"
51 #include "G4INCLNDeltaEtaProductionChannel.hh"
52 #include "G4INCLNNEtaToMultiPionsChannel.hh"
53 #include "G4INCLNNToNNOmegaChannel.hh"
54 #include "G4INCLNDeltaOmegaProductionChannel.hh"
55 #include "G4INCLNNOmegaToMultiPionsChannel.hh"
56 #include "G4INCLCrossSections.hh"
57 #include "G4INCLKinematicsUtils.hh"
58 #include "G4INCLRandom.hh"
59 #include "G4INCLParticleTable.hh"
60 #include "G4INCLPauliBlocking.hh"
61 #include "G4INCLPiNElasticChannel.hh"
62 #include "G4INCLPiNToDeltaChannel.hh"
63 #include "G4INCLPiNToMultiPionsChannel.hh"
64 #include "G4INCLPiNToEtaChannel.hh"
65 #include "G4INCLPiNToOmegaChannel.hh"
66 #include "G4INCLEtaNElasticChannel.hh"
67 #include "G4INCLEtaNToPiNChannel.hh"
68 #include "G4INCLEtaNToPiPiNChannel.hh"
69 #include "G4INCLOmegaNElasticChannel.hh"
70 #include "G4INCLOmegaNToPiNChannel.hh"
71 #include "G4INCLOmegaNToPiPiNChannel.hh"
72 #include "G4INCLStore.hh"
73 #include "G4INCLBook.hh"
74 #include "G4INCLLogger.hh"
75 #include <string>
76 #include <sstream>
77 // #include <cassert>
78 
79 namespace G4INCL {
80 
81  // WARNING: if you update the default cutNN value, make sure you update the
82  // cutNNSquared variable, too.
83  G4ThreadLocal G4double BinaryCollisionAvatar::cutNN = 1910.0;
84  G4ThreadLocal G4double BinaryCollisionAvatar::cutNNSquared = 3648100.0; // 1910.0 * 1910.0
85 
86  BinaryCollisionAvatar::BinaryCollisionAvatar(G4double time, G4double crossSection,
87  G4INCL::Nucleus *n, G4INCL::Particle *p1, G4INCL::Particle *p2)
88  : InteractionAvatar(time, n, p1, p2), theCrossSection(crossSection),
89  isParticle1Spectator(false),
90  isParticle2Spectator(false),
91  isElastic(false)
92  {
93  setType(CollisionAvatarType);
94  }
95 
96  BinaryCollisionAvatar::~BinaryCollisionAvatar() {
97  }
98 
99  G4INCL::IChannel* BinaryCollisionAvatar::getChannel() {
100  // We already check cutNN at avatar creation time, but we have to check it
101  // again here. For composite projectiles, we might have created independent
102  // avatars with no cutNN before any collision took place.
103  if(particle1->isNucleon()
104  && particle2->isNucleon()
105  && theNucleus->getStore()->getBook().getAcceptedCollisions()!=0) {
106  const G4double energyCM2 = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
107  // Below a certain cut value we don't do anything:
108  if(energyCM2 < cutNNSquared) {
109  INCL_DEBUG("CM energy = sqrt(" << energyCM2 << ") MeV < std::sqrt(" << cutNNSquared
110  << ") MeV = cutNN" << "; returning a NULL channel" << '\n');
111  InteractionAvatar::restoreParticles();
112  return NULL;
113  }
114  }
115 
135  ThreeVector minimumDistance = particle1->getPosition();
136  minimumDistance -= particle2->getPosition();
137  const G4double betaDotX = boostVector.dot(minimumDistance);
138  const G4double minDist = Math::tenPi*(minimumDistance.mag2() + betaDotX*betaDotX / (1.-boostVector.mag2()));
139  if(minDist > theCrossSection) {
140  INCL_DEBUG("CM distance of approach is too small: " << minDist << ">" <<
141  theCrossSection <<"; returning a NULL channel" << '\n');
142  InteractionAvatar::restoreParticles();
143  return NULL;
144  }
145 
147  if(particle1->isNucleon() && particle2->isNucleon()) {
148  const G4double elasticCX = CrossSections::elastic(particle1, particle2);
149  const G4double deltaProductionCX = CrossSections::NNToNDelta(particle1, particle2);
150  const G4double onePiProductionCX = CrossSections::NNToxPiNN(1,particle1, particle2);
151  const G4double twoPiProductionCX = CrossSections::NNToxPiNN(2,particle1, particle2);
152  const G4double threePiProductionCX = CrossSections::NNToxPiNN(3,particle1, particle2);
153  const G4double fourPiProductionCX = CrossSections::NNToxPiNN(4,particle1, particle2);
154  const G4double etaProductionCX = CrossSections::NNToNNEtaExclu(particle1, particle2);
155  const G4double etadeltaProductionCX = CrossSections::NNToNDeltaEta(particle1, particle2);
156  const G4double etaonePiProductionCX = CrossSections::NNToNNEtaxPi(1,particle1, particle2);
157  const G4double etatwoPiProductionCX = CrossSections::NNToNNEtaxPi(2,particle1, particle2);
158  const G4double etathreePiProductionCX = CrossSections::NNToNNEtaxPi(3,particle1, particle2);
159  const G4double etafourPiProductionCX = CrossSections::NNToNNEtaxPi(4,particle1, particle2);
160  const G4double omegaProductionCX = CrossSections::NNToNNOmegaExclu(particle1, particle2);
161  const G4double omegadeltaProductionCX = CrossSections::NNToNDeltaOmega(particle1, particle2);
162  const G4double omegaonePiProductionCX = CrossSections::NNToNNOmegaxPi(1,particle1, particle2);
163  const G4double omegatwoPiProductionCX = CrossSections::NNToNNOmegaxPi(2,particle1, particle2);
164  const G4double omegathreePiProductionCX = CrossSections::NNToNNOmegaxPi(3,particle1, particle2);
165  const G4double omegafourPiProductionCX = CrossSections::NNToNNOmegaxPi(4,particle1, particle2);
166  const G4double totCX=CrossSections::total(particle1, particle2);
167 
168  const G4double rChannel=Random::shoot() * totCX;
169 
170  if(elasticCX > rChannel) {
171 // Elastic NN channel
172  isElastic = true;
173  INCL_DEBUG("NN interaction: elastic channel chosen" << '\n');
174  return new ElasticChannel(particle1, particle2);
175  } else if((elasticCX + deltaProductionCX) > rChannel) {
176  isElastic = false;
177 // NN -> N Delta channel is chosen
178  INCL_DEBUG("NN interaction: Delta channel chosen" << '\n');
179  return new DeltaProductionChannel(particle1, particle2);
180  } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
181  isElastic = false;
182 // NN -> PiNN channel is chosen
183  INCL_DEBUG("NN interaction: one Pion channel chosen" << '\n');
184  return new NNToMultiPionsChannel(1,particle1, particle2);
185  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
186  isElastic = false;
187 // NN -> 2PiNN channel is chosen
188  INCL_DEBUG("NN interaction: two Pions channel chosen" << '\n');
189  return new NNToMultiPionsChannel(2,particle1, particle2);
190  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
191  isElastic = false;
192 // NN -> 3PiNN channel is chosen
193  INCL_DEBUG("NN interaction: three Pions channel chosen" << '\n');
194  return new NNToMultiPionsChannel(3,particle1, particle2);
195  } else if (elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX > rChannel) {
196  isElastic = false;
197 // NN -> 4PiNN channel is chosen
198  INCL_DEBUG("NN interaction: four Pions channel chosen" << '\n');
199  return new NNToMultiPionsChannel(4,particle1, particle2);
200  } else if (elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX > rChannel) {
201  isElastic = false;
202 // NN -> NNEta channel is chosen
203  INCL_DEBUG("NN interaction: Eta channel chosen" << '\n');
204  return new NNToNNEtaChannel(particle1, particle2);
205  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX > rChannel) {
206  isElastic = false;
207 // NN -> N Delta Eta channel is chosen
208  INCL_DEBUG("NN interaction: Delta Eta channel chosen" << '\n');
209  return new NDeltaEtaProductionChannel(particle1, particle2);
210  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX > rChannel) {
211  isElastic = false;
212 // NN -> EtaPiNN channel is chosen
213  INCL_DEBUG("NN interaction: Eta + one Pion channel chosen" << '\n');
214  return new NNEtaToMultiPionsChannel(1,particle1, particle2);
215  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX > rChannel) {
216  isElastic = false;
217 // NN -> Eta2PiNN channel is chosen
218  INCL_DEBUG("NN interaction: Eta + two Pions channel chosen" << '\n');
219  return new NNEtaToMultiPionsChannel(2,particle1, particle2);
220  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX > rChannel) {
221  isElastic = false;
222 // NN -> Eta3PiNN channel is chosen
223  INCL_DEBUG("NN interaction: Eta + three Pions channel chosen" << '\n');
224  return new NNEtaToMultiPionsChannel(3,particle1, particle2);
225  } else if (elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX > rChannel) {
226  isElastic = false;
227 // NN -> Eta4PiNN channel is chosen
228  INCL_DEBUG("NN interaction: Eta + four Pions channel chosen" << '\n');
229  return new NNEtaToMultiPionsChannel(4,particle1, particle2);
230  } else if (elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX + omegaProductionCX > rChannel) {
231  isElastic = false;
232 // NN -> NNOmega channel is chosen
233  INCL_DEBUG("NN interaction: Omega channel chosen" << '\n');
234  return new NNToNNOmegaChannel(particle1, particle2);
235  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX + omegaProductionCX + omegadeltaProductionCX > rChannel) {
236  isElastic = false;
237 // NN -> N Delta Omega channel is chosen
238  INCL_DEBUG("NN interaction: Delta Omega channel chosen" << '\n');
239  return new NDeltaOmegaProductionChannel(particle1, particle2);
240  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX > rChannel) {
241  isElastic = false;
242 // NN -> OmegaPiNN channel is chosen
243  INCL_DEBUG("NN interaction: Omega + one Pion channel chosen" << '\n');
244  return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
245  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX > rChannel) {
246  isElastic = false;
247 // NN -> Omega2PiNN channel is chosen
248  INCL_DEBUG("NN interaction: Omega + two Pions channel chosen" << '\n');
249  return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
250  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX > rChannel) {
251  isElastic = false;
252 // NN -> Omega3PiNN channel is chosen
253  INCL_DEBUG("NN interaction: Omega + three Pions channel chosen" << '\n');
254  return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
255  } else if (elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX > rChannel) {
256  isElastic = false;
257 // NN -> Omega4PiNN channel is chosen
258  INCL_DEBUG("NN interaction: Omega + four Pions channel chosen" << '\n');
259  return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
260  } else {
261  INCL_WARN("inconsistency within the NN Cross Sections (sum!=inelastic)" << '\n');
262  if(omegafourPiProductionCX>0.) {
263  INCL_WARN("Returning an Omega + four Pions channel" << '\n');
264  isElastic = false;
265  return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
266  } else if(omegathreePiProductionCX>0.) {
267  INCL_WARN("Returning an Omega + three Pions channel" << '\n');
268  isElastic = false;
269  return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
270  } else if(omegatwoPiProductionCX>0.) {
271  INCL_WARN("Returning an Omega + two Pions channel" << '\n');
272  isElastic = false;
273  return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
274  } else if(omegaonePiProductionCX>0.) {
275  INCL_WARN("Returning an Omega + one Pion channel" << '\n');
276  isElastic = false;
277  return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
278  } else if(omegadeltaProductionCX>0.) {
279  INCL_WARN("Returning an Omega + Delta channel" << '\n');
280  isElastic = false;
281  return new NDeltaOmegaProductionChannel(particle1, particle2);
282  } else if(omegaProductionCX>0.) {
283  INCL_WARN("Returning an Omega channel" << '\n');
284  isElastic = false;
285  return new NNToNNOmegaChannel(particle1, particle2);
286  } else if(etafourPiProductionCX>0.) {
287  INCL_WARN("Returning an Eta + four Pions channel" << '\n');
288  isElastic = false;
289  return new NNEtaToMultiPionsChannel(3,particle1, particle2);
290  } else if(etathreePiProductionCX>0.) {
291  INCL_WARN("Returning an Eta + threev channel" << '\n');
292  isElastic = false;
293  return new NNEtaToMultiPionsChannel(3,particle1, particle2);
294  } else if(etatwoPiProductionCX>0.) {
295  INCL_WARN("Returning an Eta + two Pions channel" << '\n');
296  isElastic = false;
297  return new NNEtaToMultiPionsChannel(2,particle1, particle2);
298  } else if(etaonePiProductionCX>0.) {
299  INCL_WARN("Returning an Eta + one Pion channel" << '\n');
300  isElastic = false;
301  return new NNEtaToMultiPionsChannel(1,particle1, particle2);
302  } else if(etadeltaProductionCX>0.) {
303  INCL_WARN("Returning an Eta + Delta channel" << '\n');
304  isElastic = false;
305  return new NDeltaEtaProductionChannel(particle1, particle2);
306  } else if(etaProductionCX>0.) {
307  INCL_WARN("Returning an Eta channel" << '\n');
308  isElastic = false;
309  return new NNToNNEtaChannel(particle1, particle2);
310  } else if(fourPiProductionCX>0.) {
311  INCL_WARN("Returning a 4pi channel" << '\n');
312  isElastic = false;
313  return new NNToMultiPionsChannel(4,particle1, particle2);
314  } else if(threePiProductionCX>0.) {
315  INCL_WARN("Returning a 3pi channel" << '\n');
316  isElastic = false;
317  return new NNToMultiPionsChannel(3,particle1, particle2);
318  } else if(twoPiProductionCX>0.) {
319  INCL_WARN("Returning a 2pi channel" << '\n');
320  isElastic = false;
321  return new NNToMultiPionsChannel(2,particle1, particle2);
322  } else if(onePiProductionCX>0.) {
323  INCL_WARN("Returning a 1pi channel" << '\n');
324  isElastic = false;
325  return new NNToMultiPionsChannel(1,particle1, particle2);
326  } else if(deltaProductionCX>0.) {
327  INCL_WARN("Returning a delta-production channel" << '\n');
328  isElastic = false;
329  return new DeltaProductionChannel(particle1, particle2);
330  } else {
331  INCL_WARN("Returning an elastic channel" << '\n');
332  isElastic = true;
333  return new ElasticChannel(particle1, particle2);
334  }
335  }
336 
338  } else if((particle1->isNucleon() && particle2->isDelta()) ||
339  (particle1->isDelta() && particle2->isNucleon())) {
340  G4double elasticCX = CrossSections::elastic(particle1, particle2);
341  G4double recombinationCX = CrossSections::NDeltaToNN(particle1, particle2);
342 
343  if(elasticCX/(elasticCX + recombinationCX) < Random::shoot()) {
344  isElastic = false;
345  } else
346  isElastic = true;
347 
348  if(isElastic) {
349 // Elastic N Delta channel
350  INCL_DEBUG("NDelta interaction: elastic channel chosen" << '\n');
351  return new ElasticChannel(particle1, particle2);
352  } else { // Recombination
353 // N Delta -> NN channel is chosen
354  INCL_DEBUG("NDelta interaction: recombination channel chosen" << '\n');
355  return new RecombinationChannel(particle1, particle2);
356  }
357 
359  } else if(particle1->isDelta() && particle2->isDelta()) {
360  isElastic = true;
361  INCL_DEBUG("DeltaDelta interaction: elastic channel chosen" << '\n');
362  return new ElasticChannel(particle1, particle2);
363 
365  } else if(isPiN) {
366  const G4double elasticCX = CrossSections::elastic(particle1, particle2);
367  const G4double deltaProductionCX = CrossSections::piNToDelta(particle1, particle2);
368  const G4double onePiProductionCX = CrossSections::piNToxPiN(2,particle1, particle2);
369  const G4double twoPiProductionCX = CrossSections::piNToxPiN(3,particle1, particle2);
370  const G4double threePiProductionCX = CrossSections::piNToxPiN(4,particle1, particle2);
371  const G4double etaProductionCX = CrossSections::piNToEtaN(particle1, particle2);
372  const G4double omegaProductionCX = CrossSections::piNToOmegaN(particle1, particle2);
373  const G4double totCX=CrossSections::total(particle1, particle2);
374 // assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-etaProductionCX-omegaProductionCX)<1.);
375 
376  const G4double rChannel=Random::shoot() * totCX;
377 
378  if(elasticCX > rChannel) {
379 // Elastic PiN channel
380  isElastic = true;
381  INCL_DEBUG("PiN interaction: elastic channel chosen" << '\n');
382  return new PiNElasticChannel(particle1, particle2);
383  } else if(elasticCX + deltaProductionCX > rChannel) {
384  isElastic = false;
385 // PiN -> Delta channel is chosen
386  INCL_DEBUG("PiN interaction: Delta channel chosen" << '\n');
387  return new PiNToDeltaChannel(particle1, particle2);
388  } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
389  isElastic = false;
390 // PiN -> PiNPi channel is chosen
391  INCL_DEBUG("PiN interaction: one Pion channel chosen" << '\n');
392  return new PiNToMultiPionsChannel(2,particle1, particle2);
393  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
394  isElastic = false;
395 // PiN -> PiN2Pi channel is chosen
396  INCL_DEBUG("PiN interaction: two Pions channel chosen" << '\n');
397  return new PiNToMultiPionsChannel(3,particle1, particle2);
398  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
399  isElastic = false;
400 // PiN -> PiN3Pi channel is chosen
401  INCL_DEBUG("PiN interaction: three Pions channel chosen" << '\n');
402  return new PiNToMultiPionsChannel(4,particle1, particle2);
403  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX > rChannel) {
404  isElastic = false;
405 // PiN -> EtaN channel is chosen
406  INCL_DEBUG("PiN interaction: Eta channel chosen" << '\n');
407  return new PiNToEtaChannel(particle1, particle2);
408  } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX > rChannel) {
409  isElastic = false;
410 // PiN -> OmegaN channel is chosen
411  INCL_DEBUG("PiN interaction: Omega channel chosen" << '\n');
412  return new PiNToOmegaChannel(particle1, particle2);
413  }
414 
415  else {
416  INCL_WARN("inconsistency within the PiN Cross Sections (sum!=inelastic)" << '\n');
417  if(omegaProductionCX>0.) {
418  INCL_WARN("Returning a Omega channel" << '\n');
419  isElastic = false;
420  return new PiNToOmegaChannel(particle1, particle2);
421  } else if(etaProductionCX>0.) {
422  INCL_WARN("Returning a Eta channel" << '\n');
423  isElastic = false;
424  return new PiNToEtaChannel(particle1, particle2);
425  } else if(threePiProductionCX>0.) {
426  INCL_WARN("Returning a 3pi channel" << '\n');
427  isElastic = false;
428  return new PiNToMultiPionsChannel(4,particle1, particle2);
429  } else if(twoPiProductionCX>0.) {
430  INCL_WARN("Returning a 2pi channel" << '\n');
431  isElastic = false;
432  return new PiNToMultiPionsChannel(3,particle1, particle2);
433  } else if(onePiProductionCX>0.) {
434  INCL_WARN("Returning a 1pi channel" << '\n');
435  isElastic = false;
436  return new PiNToMultiPionsChannel(2,particle1, particle2);
437  } else if(deltaProductionCX>0.) {
438  INCL_WARN("Returning a delta-production channel" << '\n');
439  isElastic = false;
440  return new PiNToDeltaChannel(particle1, particle2);
441  } else {
442  INCL_WARN("Returning an elastic channel" << '\n');
443  isElastic = true;
444  return new PiNElasticChannel(particle1, particle2);
445  }
446  }
447  } else if ((particle1->isNucleon() && particle2->isEta()) || (particle2->isNucleon() && particle1->isEta())) {
449 
450  const G4double elasticCX = CrossSections::elastic(particle1, particle2);
451  const G4double onePiProductionCX = CrossSections::etaNToPiN(particle1, particle2);
452  const G4double twoPiProductionCX = CrossSections::etaNToPiPiN(particle1, particle2);
453  const G4double totCX=CrossSections::total(particle1, particle2);
454 // assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
455 
456  const G4double rChannel=Random::shoot() * totCX;
457 
458  if(elasticCX > rChannel) {
459 // Elastic EtaN channel
460  isElastic = true;
461  INCL_DEBUG("EtaN interaction: elastic channel chosen" << '\n');
462  return new EtaNElasticChannel(particle1, particle2);
463  } else if(elasticCX + onePiProductionCX > rChannel) {
464  isElastic = false;
465 // EtaN -> EtaPiN channel is chosen
466  INCL_DEBUG("EtaN interaction: PiN channel chosen" << '\n');
467  return new EtaNToPiNChannel(particle1, particle2);
468  } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
469  isElastic = false;
470 // EtaN -> EtaPiPiN channel is chosen
471  INCL_DEBUG("EtaN interaction: PiPiN channel chosen" << '\n');
472  return new EtaNToPiPiNChannel(particle1, particle2);
473  }
474 
475  else {
476  INCL_WARN("inconsistency within the EtaN Cross Sections (sum!=inelastic)" << '\n');
477  if(twoPiProductionCX>0.) {
478  INCL_WARN("Returning a PiPiN channel" << '\n');
479  isElastic = false;
480  return new EtaNToPiPiNChannel(particle1, particle2);
481  } else if(onePiProductionCX>0.) {
482  INCL_WARN("Returning a PiN channel" << '\n');
483  isElastic = false;
484  return new EtaNToPiNChannel(particle1, particle2);
485  } else {
486  INCL_WARN("Returning an elastic channel" << '\n');
487  isElastic = true;
488  return new EtaNElasticChannel(particle1, particle2);
489  }
490  }
491 
492  } else if ((particle1->isNucleon() && particle2->isOmega()) || (particle2->isNucleon() && particle1->isOmega())) {
494 
495  const G4double elasticCX = CrossSections::elastic(particle1, particle2);
496  const G4double onePiProductionCX = CrossSections::omegaNToPiN(particle1, particle2);
497  const G4double twoPiProductionCX = CrossSections::omegaNToPiPiN(particle1, particle2);
498  const G4double totCX=CrossSections::total(particle1, particle2);
499 // assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
500 
501  const G4double rChannel=Random::shoot() * totCX;
502 
503  if(elasticCX > rChannel) {
504 // Elastic OmegaN channel
505  isElastic = true;
506  INCL_DEBUG("OmegaN interaction: elastic channel chosen" << '\n');
507  return new OmegaNElasticChannel(particle1, particle2);
508  } else if(elasticCX + onePiProductionCX > rChannel) {
509  isElastic = false;
510 // OmegaN -> OmegaPiN channel is chosen
511  INCL_DEBUG("OmegaN interaction: PiN channel chosen" << '\n');
512  return new OmegaNToPiNChannel(particle1, particle2);
513  } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
514  isElastic = false;
515 // OmegaN -> OmegaPiPiN channel is chosen
516  INCL_DEBUG("OmegaN interaction: PiPiN channel chosen" << '\n');
517  return new OmegaNToPiPiNChannel(particle1, particle2);
518  }
519 
520  else {
521  INCL_WARN("inconsistency within the OmegaN Cross Sections (sum!=inelastic)" << '\n');
522  if(twoPiProductionCX>0.) {
523  INCL_WARN("Returning a PiPiN channel" << '\n');
524  isElastic = false;
525  return new OmegaNToPiPiNChannel(particle1, particle2);
526  } else if(onePiProductionCX>0.) {
527  INCL_WARN("Returning a PiN channel" << '\n');
528  isElastic = false;
529  return new OmegaNToPiNChannel(particle1, particle2);
530  } else {
531  INCL_WARN("Returning an elastic channel" << '\n');
532  isElastic = true;
533  return new OmegaNElasticChannel(particle1, particle2);
534  }
535  }
536  }
537  else {
538  INCL_DEBUG("BinaryCollisionAvatar can only handle nucleons (for the moment)."
539  << '\n'
540  << particle1->print()
541  << '\n'
542  << particle2->print()
543  << '\n');
544  InteractionAvatar::restoreParticles();
545  return NULL;
546  }
547  }
548 
549  void BinaryCollisionAvatar::preInteraction() {
550  isParticle1Spectator = particle1->isTargetSpectator();
551  isParticle2Spectator = particle2->isTargetSpectator();
552  InteractionAvatar::preInteraction();
553  }
554 
555  void BinaryCollisionAvatar::postInteraction(FinalState *fs) {
556  // Call the postInteraction method of the parent class
557  // (provides Pauli blocking and enforces energy conservation)
558  InteractionAvatar::postInteraction(fs);
559 
560  switch(fs->getValidity()) {
561  case PauliBlockedFS:
562  theNucleus->getStore()->getBook().incrementBlockedCollisions();
563  break;
564  case NoEnergyConservationFS:
565  case ParticleBelowFermiFS:
566  case ParticleBelowZeroFS:
567  break;
568  case ValidFS:
569  Book &theBook = theNucleus->getStore()->getBook();
570  theBook.incrementAcceptedCollisions();
571  if(theBook.getAcceptedCollisions() == 1) {
572  // Store time and cross section of the first collision
573  G4double t = theBook.getCurrentTime();
574  theBook.setFirstCollisionTime(t);
575  theBook.setFirstCollisionXSec(oldXSec);
576 
577  // Store position and momentum of the spectator on the first
578  // collision
579  if((isParticle1Spectator && isParticle2Spectator) || (!isParticle1Spectator && !isParticle2Spectator)) {
580  INCL_ERROR("First collision must be within a target spectator and a non-target spectator");
581  }
582  if(isParticle1Spectator) {
583  theBook.setFirstCollisionSpectatorPosition(backupParticle1->getPosition().mag());
584  theBook.setFirstCollisionSpectatorMomentum(backupParticle1->getMomentum().mag());
585  } else {
586  theBook.setFirstCollisionSpectatorPosition(backupParticle2->getPosition().mag());
587  theBook.setFirstCollisionSpectatorMomentum(backupParticle2->getMomentum().mag());
588  }
589 
590  // Store the elasticity of the first collision
591  theBook.setFirstCollisionIsElastic(isElastic);
592  }
593  }
594  return;
595  }
596 
597  std::string BinaryCollisionAvatar::dump() const {
598  std::stringstream ss;
599  ss << "(avatar " << theTime <<" 'nn-collision" << '\n'
600  << "(list " << '\n'
601  << particle1->dump()
602  << particle2->dump()
603  << "))" << '\n';
604  return ss.str();
605  }
606 
607 }
G4INCL::Particle::isEta
G4bool isEta() const
Is this a eta?
Definition: G4INCLParticle.hh:275
G4INCLEtaNToPiPiNChannel.hh
G4INCLOmegaNElasticChannel.hh
G4INCL::CrossSections::NNToNNEtaxPi
G4double NNToNNEtaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:123
G4INCLNDeltaOmegaProductionChannel.hh
G4INCL::FinalState::getValidity
FinalStateValidity getValidity() const
Definition: G4INCLFinalState.hh:86
G4INCL::Book::incrementBlockedCollisions
void incrementBlockedCollisions()
Definition: G4INCLBook.hh:73
G4INCL::ThreeVector::dot
G4double dot(const ThreeVector &v) const
Definition: G4INCLThreeVector.hh:97
G4INCLNNEtaToMultiPionsChannel.hh
G4INCL::CrossSections::etaNToPiN
G4double etaNToPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:95
G4INCL::Book::getAcceptedCollisions
G4int getAcceptedCollisions() const
Definition: G4INCLBook.hh:100
G4INCLRandom.hh
G4INCL::KinematicsUtils::squareTotalEnergyInCM
G4double squareTotalEnergyInCM(Particle const *const p1, Particle const *const p2)
Definition: G4INCLKinematicsUtils.cc:97
G4INCL::Math::tenPi
const G4double tenPi
Definition: G4INCLGlobals.hh:70
G4INCL::CrossSections::omegaNToPiN
G4double omegaNToPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:103
G4INCLNNToNNEtaChannel.hh
INCL_ERROR
#define INCL_ERROR(x)
Definition: G4INCLLogger.hh:225
G4INCL::CrossSections::NNToNNOmegaxPi
G4double NNToNNOmegaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:139
G4INCL::Particle::isTargetSpectator
G4bool isTargetSpectator() const
Definition: G4INCLParticle.hh:251
G4INCL::CrossSections::etaNToPiPiN
G4double etaNToPiPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:99
G4INCL::Particle::getMomentum
const G4INCL::ThreeVector & getMomentum() const
Definition: G4INCLParticle.hh:578
G4INCLOmegaNToPiPiNChannel.hh
G4INCLDeltaProductionChannel.hh
G4INCL::Nucleus::getStore
Store * getStore() const
Definition: G4INCLNucleus.hh:265
G4INCL::Particle::isDelta
G4bool isDelta() const
Is it a Delta?
Definition: G4INCLParticle.hh:287
G4INCL::Particle::print
std::string print() const
Definition: G4INCLParticle.hh:737
G4INCL::CrossSections::NDeltaToNN
G4double NDeltaToNN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:63
INCL_WARN
#define INCL_WARN(x)
Definition: G4INCLLogger.hh:232
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:89
G4INCL::BinaryCollisionAvatar::postInteraction
virtual void postInteraction(FinalState *)
Definition: G4INCLBinaryCollisionAvatar.cc:555
G4INCLNNToMultiPionsChannel.hh
G4INCL::ThreeVector::mag2
G4double mag2() const
Definition: G4INCLThreeVector.hh:78
G4INCL::InteractionAvatar::backupParticle2
static G4ThreadLocal Particle * backupParticle2
Definition: G4INCLInteractionAvatar.hh:108
G4INCLBinaryCollisionAvatar.hh
G4INCL::Book::incrementAcceptedCollisions
void incrementAcceptedCollisions()
Definition: G4INCLBook.hh:72
G4INCLParticleTable.hh
G4INCL::Particle::isOmega
G4bool isOmega() const
Is this a omega?
Definition: G4INCLParticle.hh:278
G4INCL::IAvatar::setType
void setType(AvatarType t)
Definition: G4INCLIAvatar.hh:93
G4INCLPiNElasticChannel.hh
G4INCL::Store::getBook
Book & getBook()
Definition: G4INCLStore.hh:259
G4INCL::CrossSections::piNToEtaN
G4double piNToEtaN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:83
G4INCL::Book::setFirstCollisionXSec
void setFirstCollisionXSec(const G4double x)
Definition: G4INCLBook.hh:85
G4INCLPiNToMultiPionsChannel.hh
G4INCL::InteractionAvatar::backupParticle1
static G4ThreadLocal Particle * backupParticle1
Definition: G4INCLInteractionAvatar.hh:108
G4INCL::CrossSections::NNToNDeltaOmega
G4double NNToNDeltaOmega(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:143
G4INCL::CrossSections::NNToNDelta
G4double NNToNDelta(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:67
G4INCL::CrossSections::NNToNNEtaExclu
G4double NNToNNEtaExclu(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:119
G4INCLPiNToEtaChannel.hh
G4INCLOmegaNToPiNChannel.hh
G4INCL::Book::setFirstCollisionTime
void setFirstCollisionTime(const G4double t)
Definition: G4INCLBook.hh:82
G4INCLEtaNElasticChannel.hh
G4INCLNDeltaEtaProductionChannel.hh
G4INCL::CrossSections::omegaNToPiPiN
G4double omegaNToPiPiN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:107
G4INCLRecombinationChannel.hh
Delta-nucleon recombination channel.
G4INCL::CrossSections::NNToNDeltaEta
G4double NNToNDeltaEta(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:127
G4INCL::CrossSections::total
G4double total(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:59
G4INCLPiNToDeltaChannel.hh
G4INCLNNToNNOmegaChannel.hh
G4INCLCrossSections.hh
G4INCLEtaNToPiNChannel.hh
G4INCLElasticChannel.hh
G4INCLNNOmegaToMultiPionsChannel.hh
G4INCL::Particle::getPosition
const G4INCL::ThreeVector & getPosition() const
Definition: G4INCLParticle.hh:600
G4INCL::Particle::dump
std::string dump() const
Definition: G4INCLParticle.hh:752
G4INCL::CrossSections::piNToDelta
G4double piNToDelta(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:75
G4INCL::Particle::isNucleon
G4bool isNucleon() const
Definition: G4INCLParticle.hh:232
G4INCLPiNToOmegaChannel.hh
G4INCLBook.hh
G4INCL::CrossSections::NNToNNOmegaExclu
G4double NNToNNOmegaExclu(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:135
G4INCL::Random::shoot
G4double shoot()
Definition: G4INCLRandom.cc:93
G4INCL::CrossSections::piNToxPiN
G4double piNToxPiN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:79
G4INCL::InteractionAvatar::restoreParticles
void restoreParticles() const
Restore the state of both particles.
Definition: G4INCLInteractionAvatar.cc:328
G4INCL::CrossSections::piNToOmegaN
G4double piNToOmegaN(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:87
G4INCL::ThreeVector::mag
G4double mag() const
Definition: G4INCLThreeVector.hh:73
G4INCL::BinaryCollisionAvatar::BinaryCollisionAvatar
BinaryCollisionAvatar(G4double, G4double, G4INCL::Nucleus *, G4INCL::Particle *, G4INCL::Particle *)
Definition: G4INCLBinaryCollisionAvatar.cc:86
G4INCLPauliBlocking.hh
G4double
double G4double
Definition: G4Types.hh:76
G4INCL::Book::getCurrentTime
G4double getCurrentTime() const
Definition: G4INCLBook.hh:98
INCL_DEBUG
#define INCL_DEBUG(x)
Definition: G4INCLLogger.hh:240
G4INCLKinematicsUtils.hh
G4INCL::Book::setFirstCollisionIsElastic
void setFirstCollisionIsElastic(const G4bool e)
Definition: G4INCLBook.hh:94
G4INCL::CrossSections::NNToxPiNN
G4double NNToxPiNN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:71
G4INCL::Book::setFirstCollisionSpectatorMomentum
void setFirstCollisionSpectatorMomentum(const G4double x)
Definition: G4INCLBook.hh:91
G4INCLLogger.hh
G4INCL::CrossSections::elastic
G4double elastic(Particle const *const p1, Particle const *const p2)
Definition: G4INCLCrossSections.cc:55
G4INCL::Book::setFirstCollisionSpectatorPosition
void setFirstCollisionSpectatorPosition(const G4double x)
Definition: G4INCLBook.hh:88
G4INCLStore.hh