| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France 27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France
33 // 33 //
34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1
35 35
36 #include "globals.hh" 36 #include "globals.hh"
37 37
38 /* 38 /*
39 * G4INCLBinaryCollisionAvatar.cc 39 * G4INCLBinaryCollisionAvatar.cc
40 * 40 *
41 * \date Jun 5, 2009 41 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi 42 * \author Pekka Kaitaniemi
43 */ 43 */
44 44
45 #include "G4INCLBinaryCollisionAvatar.hh" 45 #include "G4INCLBinaryCollisionAvatar.hh"
46 #include "G4INCLElasticChannel.hh" 46 #include "G4INCLElasticChannel.hh"
47 #include "G4INCLRecombinationChannel.hh" 47 #include "G4INCLRecombinationChannel.hh"
48 #include "G4INCLDeltaProductionChannel.hh" 48 #include "G4INCLDeltaProductionChannel.hh"
49 #include "G4INCLNNToMultiPionsChannel.hh" 49 #include "G4INCLNNToMultiPionsChannel.hh"
50 #include "G4INCLNNToNNEtaChannel.hh" 50 #include "G4INCLNNToNNEtaChannel.hh"
51 #include "G4INCLNDeltaEtaProductionChannel.hh" 51 #include "G4INCLNDeltaEtaProductionChannel.hh"
52 #include "G4INCLNNEtaToMultiPionsChannel.hh" 52 #include "G4INCLNNEtaToMultiPionsChannel.hh"
53 #include "G4INCLNNToNNOmegaChannel.hh" 53 #include "G4INCLNNToNNOmegaChannel.hh"
54 #include "G4INCLNDeltaOmegaProductionChannel.h 54 #include "G4INCLNDeltaOmegaProductionChannel.hh"
55 #include "G4INCLNNOmegaToMultiPionsChannel.hh" 55 #include "G4INCLNNOmegaToMultiPionsChannel.hh"
56 #include "G4INCLCrossSections.hh" 56 #include "G4INCLCrossSections.hh"
57 #include "G4INCLKinematicsUtils.hh" 57 #include "G4INCLKinematicsUtils.hh"
58 #include "G4INCLRandom.hh" 58 #include "G4INCLRandom.hh"
59 #include "G4INCLParticleTable.hh" 59 #include "G4INCLParticleTable.hh"
60 #include "G4INCLPauliBlocking.hh" 60 #include "G4INCLPauliBlocking.hh"
61 #include "G4INCLPiNElasticChannel.hh" 61 #include "G4INCLPiNElasticChannel.hh"
62 #include "G4INCLPiNToDeltaChannel.hh" 62 #include "G4INCLPiNToDeltaChannel.hh"
63 #include "G4INCLPiNToMultiPionsChannel.hh" 63 #include "G4INCLPiNToMultiPionsChannel.hh"
64 #include "G4INCLPiNToEtaChannel.hh" 64 #include "G4INCLPiNToEtaChannel.hh"
65 #include "G4INCLPiNToOmegaChannel.hh" 65 #include "G4INCLPiNToOmegaChannel.hh"
66 #include "G4INCLEtaNElasticChannel.hh" 66 #include "G4INCLEtaNElasticChannel.hh"
67 #include "G4INCLEtaNToPiNChannel.hh" 67 #include "G4INCLEtaNToPiNChannel.hh"
68 #include "G4INCLEtaNToPiPiNChannel.hh" 68 #include "G4INCLEtaNToPiPiNChannel.hh"
69 #include "G4INCLOmegaNElasticChannel.hh" 69 #include "G4INCLOmegaNElasticChannel.hh"
70 #include "G4INCLOmegaNToPiNChannel.hh" 70 #include "G4INCLOmegaNToPiNChannel.hh"
71 #include "G4INCLNNToNLKChannel.hh" 71 #include "G4INCLNNToNLKChannel.hh"
72 #include "G4INCLNNToNSKChannel.hh" 72 #include "G4INCLNNToNSKChannel.hh"
73 #include "G4INCLNNToNLKpiChannel.hh" 73 #include "G4INCLNNToNLKpiChannel.hh"
74 #include "G4INCLNNToNSKpiChannel.hh" 74 #include "G4INCLNNToNSKpiChannel.hh"
75 #include "G4INCLNNToNLK2piChannel.hh" 75 #include "G4INCLNNToNLK2piChannel.hh"
76 #include "G4INCLNNToNSK2piChannel.hh" 76 #include "G4INCLNNToNSK2piChannel.hh"
77 #include "G4INCLNNToNNKKbChannel.hh" 77 #include "G4INCLNNToNNKKbChannel.hh"
78 #include "G4INCLNNToMissingStrangenessChannel. 78 #include "G4INCLNNToMissingStrangenessChannel.hh"
79 #include "G4INCLNDeltaToNLKChannel.hh" 79 #include "G4INCLNDeltaToNLKChannel.hh"
80 #include "G4INCLNDeltaToNSKChannel.hh" 80 #include "G4INCLNDeltaToNSKChannel.hh"
81 #include "G4INCLNDeltaToDeltaLKChannel.hh" 81 #include "G4INCLNDeltaToDeltaLKChannel.hh"
82 #include "G4INCLNDeltaToDeltaSKChannel.hh" 82 #include "G4INCLNDeltaToDeltaSKChannel.hh"
83 #include "G4INCLNDeltaToNNKKbChannel.hh" 83 #include "G4INCLNDeltaToNNKKbChannel.hh"
84 #include "G4INCLNpiToLKChannel.hh" 84 #include "G4INCLNpiToLKChannel.hh"
85 #include "G4INCLNpiToSKChannel.hh" 85 #include "G4INCLNpiToSKChannel.hh"
86 #include "G4INCLNpiToLKpiChannel.hh" 86 #include "G4INCLNpiToLKpiChannel.hh"
87 #include "G4INCLNpiToSKpiChannel.hh" 87 #include "G4INCLNpiToSKpiChannel.hh"
88 #include "G4INCLNpiToLK2piChannel.hh" 88 #include "G4INCLNpiToLK2piChannel.hh"
89 #include "G4INCLNpiToSK2piChannel.hh" 89 #include "G4INCLNpiToSK2piChannel.hh"
90 #include "G4INCLNpiToNKKbChannel.hh" 90 #include "G4INCLNpiToNKKbChannel.hh"
91 #include "G4INCLNpiToMissingStrangenessChannel 91 #include "G4INCLNpiToMissingStrangenessChannel.hh"
92 #include "G4INCLNKElasticChannel.hh" 92 #include "G4INCLNKElasticChannel.hh"
93 #include "G4INCLNKToNKChannel.hh" 93 #include "G4INCLNKToNKChannel.hh"
94 #include "G4INCLNKToNKpiChannel.hh" 94 #include "G4INCLNKToNKpiChannel.hh"
95 #include "G4INCLNKToNK2piChannel.hh" 95 #include "G4INCLNKToNK2piChannel.hh"
96 #include "G4INCLNKbElasticChannel.hh" 96 #include "G4INCLNKbElasticChannel.hh"
97 #include "G4INCLNKbToNKbChannel.hh" 97 #include "G4INCLNKbToNKbChannel.hh"
98 #include "G4INCLNKbToNKbpiChannel.hh" 98 #include "G4INCLNKbToNKbpiChannel.hh"
99 #include "G4INCLNKbToNKb2piChannel.hh" 99 #include "G4INCLNKbToNKb2piChannel.hh"
100 #include "G4INCLNKbToLpiChannel.hh" 100 #include "G4INCLNKbToLpiChannel.hh"
101 #include "G4INCLNKbToL2piChannel.hh" 101 #include "G4INCLNKbToL2piChannel.hh"
102 #include "G4INCLNKbToSpiChannel.hh" 102 #include "G4INCLNKbToSpiChannel.hh"
103 #include "G4INCLNKbToS2piChannel.hh" 103 #include "G4INCLNKbToS2piChannel.hh"
104 #include "G4INCLNYElasticChannel.hh" 104 #include "G4INCLNYElasticChannel.hh"
105 #include "G4INCLNLToNSChannel.hh" 105 #include "G4INCLNLToNSChannel.hh"
106 #include "G4INCLNSToNLChannel.hh" 106 #include "G4INCLNSToNLChannel.hh"
107 #include "G4INCLNSToNSChannel.hh" 107 #include "G4INCLNSToNSChannel.hh"
108 #include "G4INCLOmegaNToPiPiNChannel.hh" 108 #include "G4INCLOmegaNToPiPiNChannel.hh"
109 #include "G4INCLStore.hh" 109 #include "G4INCLStore.hh"
110 #include "G4INCLBook.hh" 110 #include "G4INCLBook.hh"
111 #include "G4INCLLogger.hh" 111 #include "G4INCLLogger.hh"
112 #include <string> 112 #include <string>
113 #include <sstream> 113 #include <sstream>
114 // #include <cassert> 114 // #include <cassert>
115 #include "G4INCLNNbarElasticChannel.hh" <<
116 #include "G4INCLNNbarCEXChannel.hh" <<
117 #include "G4INCLNNbarToLLbarChannel.hh" <<
118 #include "G4INCLNNbarToNNbarpiChannel.hh" <<
119 #include "G4INCLNNbarToNNbar2piChannel.hh" <<
120 #include "G4INCLNNbarToNNbar3piChannel.hh" <<
121 #include "G4INCLNNbarToAnnihilationChannel.hh" <<
122 115
123 namespace G4INCL { 116 namespace G4INCL {
124 117
125 // WARNING: if you update the default cutNN 118 // WARNING: if you update the default cutNN value, make sure you update the
126 // cutNNSquared variable, too. 119 // cutNNSquared variable, too.
127 G4ThreadLocal G4double BinaryCollisionAvatar 120 G4ThreadLocal G4double BinaryCollisionAvatar::cutNN = 1910.0;
128 G4ThreadLocal G4double BinaryCollisionAvatar 121 G4ThreadLocal G4double BinaryCollisionAvatar::cutNNSquared = 3648100.0; // 1910.0 * 1910.0
129 G4ThreadLocal G4double BinaryCollisionAvatar 122 G4ThreadLocal G4double BinaryCollisionAvatar::bias = 1.;
130 123
131 BinaryCollisionAvatar::BinaryCollisionAvatar 124 BinaryCollisionAvatar::BinaryCollisionAvatar(G4double time, G4double crossSection,
132 G4INCL::Nucleus *n, G4INCL::Particle *p1 125 G4INCL::Nucleus *n, G4INCL::Particle *p1, G4INCL::Particle *p2)
133 : InteractionAvatar(time, n, p1, p2), theC 126 : InteractionAvatar(time, n, p1, p2), theCrossSection(crossSection),
134 isParticle1Spectator(false), 127 isParticle1Spectator(false),
135 isParticle2Spectator(false), 128 isParticle2Spectator(false),
136 isElastic(false), 129 isElastic(false),
137 isStrangeProduction(false) 130 isStrangeProduction(false)
138 { 131 {
139 setType(CollisionAvatarType); 132 setType(CollisionAvatarType);
140 } 133 }
141 134
142 BinaryCollisionAvatar::~BinaryCollisionAvata 135 BinaryCollisionAvatar::~BinaryCollisionAvatar() {
143 } 136 }
144 137
145 G4INCL::IChannel* BinaryCollisionAvatar::get 138 G4INCL::IChannel* BinaryCollisionAvatar::getChannel() {
146 // We already check cutNN at avatar creati 139 // We already check cutNN at avatar creation time, but we have to check it
147 // again here. For composite projectiles, 140 // again here. For composite projectiles, we might have created independent
148 // avatars with no cutNN before any collis 141 // avatars with no cutNN before any collision took place.
149 if(particle1->isNucleon() 142 if(particle1->isNucleon()
150 && particle2->isNucleon() 143 && particle2->isNucleon()
151 && theNucleus->getStore()->getBook().g 144 && theNucleus->getStore()->getBook().getAcceptedCollisions()!=0) {
152 const G4double energyCM2 = KinematicsUti 145 const G4double energyCM2 = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2);
153 // Below a certain cut value we don't do 146 // Below a certain cut value we don't do anything:
154 if(energyCM2 < cutNNSquared) { 147 if(energyCM2 < cutNNSquared) {
155 INCL_DEBUG("CM energy = sqrt(" << ener 148 INCL_DEBUG("CM energy = sqrt(" << energyCM2 << ") MeV < std::sqrt(" << cutNNSquared
156 << ") MeV = cutNN" << "; returning 149 << ") MeV = cutNN" << "; returning a NULL channel" << '\n');
157 InteractionAvatar::restoreParticles(); 150 InteractionAvatar::restoreParticles();
158 return NULL; 151 return NULL;
159 } 152 }
160 } 153 }
161 154
162 /** Check again the distance of approach. 155 /** Check again the distance of approach. In order for the avatar to be
163 * realised, we have to perform a check in 156 * realised, we have to perform a check in the CM system. We define a
164 * distance four-vector as 157 * distance four-vector as
165 * \f[ (0, \Delta\vec{x}), \f] 158 * \f[ (0, \Delta\vec{x}), \f]
166 * where \f$\Delta\vec{x}\f$ is the distan 159 * where \f$\Delta\vec{x}\f$ is the distance vector of the particles at
167 * their minimum distance of approach (i.e 160 * their minimum distance of approach (i.e. at the avatar time). By
168 * boosting this four-vector to the CM fra 161 * boosting this four-vector to the CM frame of the two particles and we
169 * obtain a new four vector 162 * obtain a new four vector
170 * \f[ (\Delta t', \Delta\vec{x}'), \f] 163 * \f[ (\Delta t', \Delta\vec{x}'), \f]
171 * with a non-zero time component (the col 164 * with a non-zero time component (the collision happens simultaneously for
172 * the two particles in the lab system, bu 165 * the two particles in the lab system, but not in the CM system). In order
173 * for the avatar to be realised, we requi 166 * for the avatar to be realised, we require that
174 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/ 167 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/\pi}.\f]
175 * Note that \f$|\Delta\vec{x}'|\leq|\Delt 168 * Note that \f$|\Delta\vec{x}'|\leq|\Delta\vec{x}|\f$; thus, the condition
176 * above is more restrictive than the chec 169 * above is more restrictive than the check that we perform in
177 * G4INCL::Propagation::StandardPropagatio 170 * G4INCL::Propagation::StandardPropagationModel::generateBinaryCollisionAvatar.
178 * In other words, the avatar generation c 171 * In other words, the avatar generation cannot miss any physical collision
179 * avatars. 172 * avatars.
180 */ 173 */
181 ThreeVector minimumDistance = particle1->g 174 ThreeVector minimumDistance = particle1->getPosition();
182 minimumDistance -= particle2->getPosition( 175 minimumDistance -= particle2->getPosition();
183 const G4double betaDotX = boostVector.dot( 176 const G4double betaDotX = boostVector.dot(minimumDistance);
184 const G4double minDist = Math::tenPi*(mini 177 const G4double minDist = Math::tenPi*(minimumDistance.mag2() + betaDotX*betaDotX / (1.-boostVector.mag2()));
185 if(minDist > theCrossSection) { 178 if(minDist > theCrossSection) {
186 INCL_DEBUG("CM distance of approach is t 179 INCL_DEBUG("CM distance of approach is too small: " << minDist << ">" <<
187 theCrossSection <<"; returning a NULL 180 theCrossSection <<"; returning a NULL channel" << '\n');
188 InteractionAvatar::restoreParticles(); 181 InteractionAvatar::restoreParticles();
189 return NULL; 182 return NULL;
190 } 183 }
191 184
192 /** Bias apply for this reaction in order 185 /** Bias apply for this reaction in order to get the same
193 * ParticleBias for all stange particles. 186 * ParticleBias for all stange particles.
194 * Can be reduced after because of the saf << 187 * Can be reduced after because of the watchdog.
195 */ 188 */
196 G4double bias_apply = 1.; 189 G4double bias_apply = 1.;
197 if(bias != 1.) bias_apply = Particle::getB 190 if(bias != 1.) bias_apply = Particle::getBiasFromVector(Particle::MergeVectorBias(particle1,particle2)) * bias;
198 191
199 //// NN 192 //// NN
200 if(particle1->isNucleon() && particle2->is 193 if(particle1->isNucleon() && particle2->isNucleon()) {
201 194
202 G4double NLKProductionCX = CrossSections 195 G4double NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
203 G4double NSKProductionCX = CrossSections 196 G4double NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
204 G4double NLKpiProductionCX = CrossSectio 197 G4double NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
205 G4double NSKpiProductionCX = CrossSectio 198 G4double NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
206 G4double NLK2piProductionCX = CrossSecti 199 G4double NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
207 G4double NSK2piProductionCX = CrossSecti 200 G4double NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
208 G4double NNKKbProductionCX = CrossSectio 201 G4double NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
209 G4double NNMissingCX = CrossSections::NN 202 G4double NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
210 203
211 const G4double UnStrangeProdCX = CrossSe 204 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NNToNDelta(particle1, particle2) + CrossSections::NNToxPiNN(1,particle1, particle2)
212 + CrossSect 205 + CrossSections::NNToxPiNN(2,particle1, particle2) + CrossSections::NNToxPiNN(3,particle1, particle2) + CrossSections::NNToxPiNN(4,particle1, particle2)
213 + CrossSect 206 + CrossSections::NNToNNEtaExclu(particle1, particle2) + CrossSections::NNToNDeltaEta(particle1, particle2) + CrossSections::NNToNNEtaxPi(1,particle1, particle2)
214 + CrossSect 207 + CrossSections::NNToNNEtaxPi(2,particle1, particle2) + CrossSections::NNToNNEtaxPi(3,particle1, particle2) + CrossSections::NNToNNEtaxPi(4,particle1, particle2)
215 + CrossSect 208 + CrossSections::NNToNNOmegaExclu(particle1, particle2) + CrossSections::NNToNDeltaOmega(particle1, particle2) + CrossSections::NNToNNOmegaxPi(1,particle1, particle2)
216 + CrossSect 209 + CrossSections::NNToNNOmegaxPi(2,particle1, particle2) + CrossSections::NNToNNOmegaxPi(3,particle1, particle2) + CrossSections::NNToNNOmegaxPi(4,particle1, particle2);
217 const G4double StrangenessProdCX = (NLKP 210 const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + NLKpiProductionCX + NSKpiProductionCX + NLK2piProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX)/bias_apply;
218 211
219 G4double counterweight = (1. - bias_appl 212 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
220 213
221 if(counterweight < 0.5) { 214 if(counterweight < 0.5) {
222 counterweight = 0.5; 215 counterweight = 0.5;
223 bias_apply = 0.5*UnStrangeProdCX/Stra 216 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
224 NLKProductionCX = CrossSections::NNTo 217 NLKProductionCX = CrossSections::NNToNLK(particle1, particle2)*bias_apply;
225 NSKProductionCX = CrossSections::NNTo 218 NSKProductionCX = CrossSections::NNToNSK(particle1, particle2)*bias_apply;
226 NLKpiProductionCX = CrossSections::NN 219 NLKpiProductionCX = CrossSections::NNToNLKpi(particle1, particle2)*bias_apply;
227 NSKpiProductionCX = CrossSections::NN 220 NSKpiProductionCX = CrossSections::NNToNSKpi(particle1, particle2)*bias_apply;
228 NLK2piProductionCX = CrossSections::N 221 NLK2piProductionCX = CrossSections::NNToNLK2pi(particle1, particle2)*bias_apply;
229 NSK2piProductionCX = CrossSections::N 222 NSK2piProductionCX = CrossSections::NNToNSK2pi(particle1, particle2)*bias_apply;
230 NNKKbProductionCX = CrossSections::NN 223 NNKKbProductionCX = CrossSections::NNToNNKKb(particle1, particle2)*bias_apply;
231 NNMissingCX = CrossSections::NNToMiss 224 NNMissingCX = CrossSections::NNToMissingStrangeness(particle1, particle2)*bias_apply;
232 } 225 }
233 226
234 227
235 const G4double elasticCX = CrossSections 228 const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
236 const G4double deltaProductionCX = Cross 229 const G4double deltaProductionCX = CrossSections::NNToNDelta(particle1, particle2)*counterweight;
237 const G4double onePiProductionCX = Cross 230 const G4double onePiProductionCX = CrossSections::NNToxPiNN(1,particle1, particle2)*counterweight;
238 const G4double twoPiProductionCX = Cross 231 const G4double twoPiProductionCX = CrossSections::NNToxPiNN(2,particle1, particle2)*counterweight;
239 const G4double threePiProductionCX = Cro 232 const G4double threePiProductionCX = CrossSections::NNToxPiNN(3,particle1, particle2)*counterweight;
240 const G4double fourPiProductionCX = Cros 233 const G4double fourPiProductionCX = CrossSections::NNToxPiNN(4,particle1, particle2)*counterweight;
241 234
242 const G4double etaProductionCX = CrossSe 235 const G4double etaProductionCX = CrossSections::NNToNNEtaExclu(particle1, particle2)*counterweight;
243 const G4double etadeltaProductionCX = Cr 236 const G4double etadeltaProductionCX = CrossSections::NNToNDeltaEta(particle1, particle2)*counterweight;
244 const G4double etaonePiProductionCX = Cr 237 const G4double etaonePiProductionCX = CrossSections::NNToNNEtaxPi(1,particle1, particle2)*counterweight;
245 const G4double etatwoPiProductionCX = Cr 238 const G4double etatwoPiProductionCX = CrossSections::NNToNNEtaxPi(2,particle1, particle2)*counterweight;
246 const G4double etathreePiProductionCX = 239 const G4double etathreePiProductionCX = CrossSections::NNToNNEtaxPi(3,particle1, particle2)*counterweight;
247 const G4double etafourPiProductionCX = C 240 const G4double etafourPiProductionCX = CrossSections::NNToNNEtaxPi(4,particle1, particle2)*counterweight;
248 const G4double omegaProductionCX = Cross 241 const G4double omegaProductionCX = CrossSections::NNToNNOmegaExclu(particle1, particle2)*counterweight;
249 const G4double omegadeltaProductionCX = 242 const G4double omegadeltaProductionCX = CrossSections::NNToNDeltaOmega(particle1, particle2)*counterweight;
250 const G4double omegaonePiProductionCX = 243 const G4double omegaonePiProductionCX = CrossSections::NNToNNOmegaxPi(1,particle1, particle2)*counterweight;
251 const G4double omegatwoPiProductionCX = 244 const G4double omegatwoPiProductionCX = CrossSections::NNToNNOmegaxPi(2,particle1, particle2)*counterweight;
252 const G4double omegathreePiProductionCX 245 const G4double omegathreePiProductionCX = CrossSections::NNToNNOmegaxPi(3,particle1, particle2)*counterweight;
253 const G4double omegafourPiProductionCX = 246 const G4double omegafourPiProductionCX = CrossSections::NNToNNOmegaxPi(4,particle1, particle2)*counterweight;
254 247
255 const G4double totCX=CrossSections::tota 248 const G4double totCX=CrossSections::total(particle1, particle2);
256 249
257 // assert(std::fabs(totCX-elasticCX-deltaProdu 250 // assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-fourPiProductionCX-NLKProductionCX-NSKProductionCX-NLKpiProductionCX-NSKpiProductionCX-NLK2piProductionCX-NSK2piProductionCX-NNKKbProductionCX-NNMissingCX-etaProductionCX-etadeltaProductionCX-etaonePiProductionCX-etatwoPiProductionCX-etathreePiProductionCX-etafourPiProductionCX-omegaProductionCX-omegadeltaProductionCX-omegaonePiProductionCX-omegatwoPiProductionCX-omegathreePiProductionCX-omegafourPiProductionCX) < 0.5);
258 251
259 const G4double rChannel=Random::shoot() 252 const G4double rChannel=Random::shoot() * totCX;
260 253
261 if(elasticCX > rChannel) { 254 if(elasticCX > rChannel) {
262 // Elastic NN channel 255 // Elastic NN channel
263 isElastic = true; 256 isElastic = true;
264 INCL_DEBUG("NN interaction: elastic ch 257 INCL_DEBUG("NN interaction: elastic channel chosen" << '\n');
265 weight = counterweight; 258 weight = counterweight;
266 return new ElasticChannel(particle1, p 259 return new ElasticChannel(particle1, particle2);
267 } else if((elasticCX + deltaProductionCX 260 } else if((elasticCX + deltaProductionCX) > rChannel) {
268 isElastic = false; 261 isElastic = false;
269 // NN -> N Delta channel is chosen 262 // NN -> N Delta channel is chosen
270 INCL_DEBUG("NN interaction: Delta chan 263 INCL_DEBUG("NN interaction: Delta channel chosen" << '\n');
271 weight = counterweight; 264 weight = counterweight;
272 return new DeltaProductionChannel(part 265 return new DeltaProductionChannel(particle1, particle2);
273 } else if(elasticCX + deltaProductionCX 266 } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
274 isElastic = false; 267 isElastic = false;
275 // NN -> PiNN channel is chosen 268 // NN -> PiNN channel is chosen
276 INCL_DEBUG("NN interaction: one Pion c 269 INCL_DEBUG("NN interaction: one Pion channel chosen" << '\n');
277 weight = counterweight; 270 weight = counterweight;
278 return new NNToMultiPionsChannel(1,par 271 return new NNToMultiPionsChannel(1,particle1, particle2);
279 } else if(elasticCX + deltaProductionCX 272 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
280 isElastic = false; 273 isElastic = false;
281 // NN -> 2PiNN channel is chosen 274 // NN -> 2PiNN channel is chosen
282 INCL_DEBUG("NN interaction: two Pions 275 INCL_DEBUG("NN interaction: two Pions channel chosen" << '\n');
283 weight = counterweight; 276 weight = counterweight;
284 return new NNToMultiPionsChannel(2,par 277 return new NNToMultiPionsChannel(2,particle1, particle2);
285 } else if(elasticCX + deltaProductionCX 278 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
286 isElastic = false; 279 isElastic = false;
287 // NN -> 3PiNN channel is chosen 280 // NN -> 3PiNN channel is chosen
288 INCL_DEBUG("NN interaction: three Pion 281 INCL_DEBUG("NN interaction: three Pions channel chosen" << '\n');
289 weight = counterweight; 282 weight = counterweight;
290 return new NNToMultiPionsChannel(3,par 283 return new NNToMultiPionsChannel(3,particle1, particle2);
291 } else if(elasticCX + deltaProductionCX 284 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX > rChannel) {
292 isElastic = false; 285 isElastic = false;
293 // NN -> 4PiNN channel is chosen 286 // NN -> 4PiNN channel is chosen
294 INCL_DEBUG("NN interaction: four 287 INCL_DEBUG("NN interaction: four Pions channel chosen" << '\n');
295 weight = counterweight; 288 weight = counterweight;
296 return new NNToMultiPionsChannel 289 return new NNToMultiPionsChannel(4,particle1, particle2);
297 } else if(elasticCX + deltaProductionCX 290 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
298 + etaProductionCX > 291 + etaProductionCX > rChannel) {
299 isElastic = false; 292 isElastic = false;
300 // NN -> NNEta channel is chosen 293 // NN -> NNEta channel is chosen
301 INCL_DEBUG("NN interaction: Eta channel 294 INCL_DEBUG("NN interaction: Eta channel chosen" << '\n');
302 weight = counterweight; 295 weight = counterweight;
303 return new NNToNNEtaChannel(particle1, 296 return new NNToNNEtaChannel(particle1, particle2);
304 } else if(elasticCX + deltaProductionCX 297 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
305 + etaProductionCX + 298 + etaProductionCX + etadeltaProductionCX > rChannel) {
306 isElastic = false; 299 isElastic = false;
307 // NN -> N Delta Eta channel is chosen 300 // NN -> N Delta Eta channel is chosen
308 INCL_DEBUG("NN interaction: Delta Eta c 301 INCL_DEBUG("NN interaction: Delta Eta channel chosen" << '\n');
309 weight = counterweight; 302 weight = counterweight;
310 return new NDeltaEtaProductionChannel(p 303 return new NDeltaEtaProductionChannel(particle1, particle2);
311 } else if(elasticCX + deltaProductionCX 304 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
312 + etaProductionCX + 305 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX > rChannel) {
313 isElastic = false; 306 isElastic = false;
314 // NN -> EtaPiNN channel is chosen 307 // NN -> EtaPiNN channel is chosen
315 INCL_DEBUG("NN interaction: Eta + one P 308 INCL_DEBUG("NN interaction: Eta + one Pion channel chosen" << '\n');
316 weight = counterweight; 309 weight = counterweight;
317 return new NNEtaToMultiPionsChannel(1,p 310 return new NNEtaToMultiPionsChannel(1,particle1, particle2);
318 } else if(elasticCX + deltaProductionCX 311 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
319 + etaProductionCX + 312 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX > rChannel) {
320 isElastic = false; 313 isElastic = false;
321 // NN -> Eta2PiNN channel is chosen 314 // NN -> Eta2PiNN channel is chosen
322 INCL_DEBUG("NN interaction: Eta + two P 315 INCL_DEBUG("NN interaction: Eta + two Pions channel chosen" << '\n');
323 weight = counterweight; 316 weight = counterweight;
324 return new NNEtaToMultiPionsChannel(2,p 317 return new NNEtaToMultiPionsChannel(2,particle1, particle2);
325 } else if(elasticCX + deltaProductionCX 318 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
326 + etaProductionCX + 319 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX > rChannel) {
327 isElastic = false; 320 isElastic = false;
328 // NN -> Eta3PiNN channel is chosen 321 // NN -> Eta3PiNN channel is chosen
329 INCL_DEBUG("NN interaction: Eta + three 322 INCL_DEBUG("NN interaction: Eta + three Pions channel chosen" << '\n');
330 weight = counterweight; 323 weight = counterweight;
331 return new NNEtaToMultiPionsChannel(3,p 324 return new NNEtaToMultiPionsChannel(3,particle1, particle2);
332 } else if(elasticCX + deltaProductionCX 325 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
333 + etaProductionCX + 326 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX > rChannel) {
334 isElastic = false; 327 isElastic = false;
335 // NN -> Eta4PiNN channel is chosen 328 // NN -> Eta4PiNN channel is chosen
336 INCL_DEBUG("NN interaction: Eta + four 329 INCL_DEBUG("NN interaction: Eta + four Pions channel chosen" << '\n');
337 weight = counterweight; 330 weight = counterweight;
338 return new NNEtaToMultiPionsChannel(4,p 331 return new NNEtaToMultiPionsChannel(4,particle1, particle2);
339 } else if(elasticCX + deltaProductionCX 332 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
340 + etaProductionCX + 333 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
341 + omegaProductionCX 334 + omegaProductionCX > rChannel) {
342 isElastic = false; 335 isElastic = false;
343 // NN -> NNOmega channel is chosen 336 // NN -> NNOmega channel is chosen
344 INCL_DEBUG("NN interaction: Omega chann 337 INCL_DEBUG("NN interaction: Omega channel chosen" << '\n');
345 weight = counterweight; 338 weight = counterweight;
346 return new NNToNNOmegaChannel(particle1 339 return new NNToNNOmegaChannel(particle1, particle2);
347 } else if(elasticCX + deltaProductionCX 340 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
348 + etaProductionCX + 341 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
349 + omegaProductionCX 342 + omegaProductionCX + omegadeltaProductionCX > rChannel) {
350 isElastic = false; 343 isElastic = false;
351 // NN -> N Delta Omega channel is chosen 344 // NN -> N Delta Omega channel is chosen
352 INCL_DEBUG("NN interaction: Delta Omega 345 INCL_DEBUG("NN interaction: Delta Omega channel chosen" << '\n');
353 weight = counterweight; 346 weight = counterweight;
354 return new NDeltaOmegaProductionChannel 347 return new NDeltaOmegaProductionChannel(particle1, particle2);
355 } else if(elasticCX + deltaProductionCX 348 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
356 + etaProductionCX + 349 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
357 + omegaProductionCX 350 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX > rChannel) {
358 isElastic = false; 351 isElastic = false;
359 // NN -> OmegaPiNN channel is chosen 352 // NN -> OmegaPiNN channel is chosen
360 INCL_DEBUG("NN interaction: Omega + one 353 INCL_DEBUG("NN interaction: Omega + one Pion channel chosen" << '\n');
361 weight = counterweight; 354 weight = counterweight;
362 return new NNOmegaToMultiPionsChannel(1 355 return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
363 } else if(elasticCX + deltaProductionCX 356 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
364 + etaProductionCX + 357 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
365 + omegaProductionCX 358 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX > rChannel) {
366 isElastic = false; 359 isElastic = false;
367 // NN -> Omega2PiNN channel is chosen 360 // NN -> Omega2PiNN channel is chosen
368 INCL_DEBUG("NN interaction: Omega + two 361 INCL_DEBUG("NN interaction: Omega + two Pions channel chosen" << '\n');
369 weight = counterweight; 362 weight = counterweight;
370 return new NNOmegaToMultiPionsChannel(2 363 return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
371 } else if(elasticCX + deltaProductionCX 364 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
372 + etaProductionCX + 365 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
373 + omegaProductionCX 366 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX > rChannel) {
374 isElastic = false; 367 isElastic = false;
375 // NN -> Omega3PiNN channel is chosen 368 // NN -> Omega3PiNN channel is chosen
376 INCL_DEBUG("NN interaction: Omega + thr 369 INCL_DEBUG("NN interaction: Omega + three Pions channel chosen" << '\n');
377 weight = counterweight; 370 weight = counterweight;
378 return new NNOmegaToMultiPionsChannel(3 371 return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
379 } else if(elasticCX + deltaProductionCX 372 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
380 + etaProductionCX + 373 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
381 + omegaProductionCX 374 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX > rChannel) {
382 isElastic = false; 375 isElastic = false;
383 // NN -> Omega4PiNN channel is chosen 376 // NN -> Omega4PiNN channel is chosen
384 INCL_DEBUG("NN interaction: Omega + fou 377 INCL_DEBUG("NN interaction: Omega + four Pions channel chosen" << '\n');
385 weight = counterweight; 378 weight = counterweight;
386 return new NNOmegaToMultiPionsChannel(4 379 return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
387 } else if(elasticCX + deltaProductionCX 380 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
388 + etaProductionCX + 381 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
389 + omegaProductionCX 382 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
390 + NLKProductionCX > 383 + NLKProductionCX > rChannel) {
391 isElastic = false; 384 isElastic = false;
392 isStrangeProduction = true; 385 isStrangeProduction = true;
393 // NN -> NLK channel is chosen 386 // NN -> NLK channel is chosen
394 INCL_DEBUG("NN interaction: NLK channe 387 INCL_DEBUG("NN interaction: NLK channel chosen" << '\n');
395 weight = bias_apply; 388 weight = bias_apply;
396 return new NNToNLKChannel(particle1, p 389 return new NNToNLKChannel(particle1, particle2);
397 } else if(elasticCX + deltaProductionCX 390 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
398 + etaProductionCX + 391 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
399 + omegaProductionCX 392 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
400 + NLKProductionCX + 393 + NLKProductionCX + NLKpiProductionCX > rChannel) {
401 isElastic = false; 394 isElastic = false;
402 isStrangeProduction = true; 395 isStrangeProduction = true;
403 // NN -> NLKpi channel is chosen 396 // NN -> NLKpi channel is chosen
404 INCL_DEBUG("NN interaction: NLKpi chan 397 INCL_DEBUG("NN interaction: NLKpi channel chosen" << '\n');
405 weight = bias_apply; 398 weight = bias_apply;
406 return new NNToNLKpiChannel(particle1, 399 return new NNToNLKpiChannel(particle1, particle2);
407 } else if(elasticCX + deltaProductionCX 400 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
408 + etaProductionCX + 401 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
409 + omegaProductionCX 402 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
410 + NLKProductionCX + 403 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX > rChannel) {
411 isElastic = false; 404 isElastic = false;
412 isStrangeProduction = true; 405 isStrangeProduction = true;
413 // NN -> NLK2pi channel is chosen 406 // NN -> NLK2pi channel is chosen
414 INCL_DEBUG("NN interaction: NLK2pi cha 407 INCL_DEBUG("NN interaction: NLK2pi channel chosen" << '\n');
415 weight = bias_apply; 408 weight = bias_apply;
416 return new NNToNLK2piChannel(particle1 409 return new NNToNLK2piChannel(particle1, particle2);
417 } else if(elasticCX + deltaProductionCX 410 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
418 + etaProductionCX + 411 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
419 + omegaProductionCX 412 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
420 + NLKProductionCX + 413 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX > rChannel) {
421 isElastic = false; 414 isElastic = false;
422 isStrangeProduction = true; 415 isStrangeProduction = true;
423 // NN -> NSK channel is chosen 416 // NN -> NSK channel is chosen
424 INCL_DEBUG("NN interaction: NSK channe 417 INCL_DEBUG("NN interaction: NSK channel chosen" << '\n');
425 weight = bias_apply; 418 weight = bias_apply;
426 return new NNToNSKChannel(particle1, p 419 return new NNToNSKChannel(particle1, particle2);
427 } else if(elasticCX + deltaProductionCX 420 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
428 + etaProductionCX + 421 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
429 + omegaProductionCX 422 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
430 + NLKProductionCX + 423 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX > rChannel) {
431 isElastic = false; 424 isElastic = false;
432 isStrangeProduction = true; 425 isStrangeProduction = true;
433 // NN -> NSKpi channel is chosen 426 // NN -> NSKpi channel is chosen
434 INCL_DEBUG("NN interaction: NSKpi chan 427 INCL_DEBUG("NN interaction: NSKpi channel chosen" << '\n');
435 weight = bias_apply; 428 weight = bias_apply;
436 return new NNToNSKpiChannel(particle1, 429 return new NNToNSKpiChannel(particle1, particle2);
437 } else if(elasticCX + deltaProductionCX 430 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
438 + etaProductionCX + 431 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
439 + omegaProductionCX 432 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
440 + NLKProductionCX + 433 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX > rChannel) {
441 isElastic = false; 434 isElastic = false;
442 isStrangeProduction = true; 435 isStrangeProduction = true;
443 // NN -> NSK2pi channel is chosen 436 // NN -> NSK2pi channel is chosen
444 INCL_DEBUG("NN interaction: NSK2pi cha 437 INCL_DEBUG("NN interaction: NSK2pi channel chosen" << '\n');
445 weight = bias_apply; 438 weight = bias_apply;
446 return new NNToNSK2piChannel(particle1 439 return new NNToNSK2piChannel(particle1, particle2);
447 } else if(elasticCX + deltaProductionCX 440 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
448 + etaProductionCX + 441 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
449 + omegaProductionCX 442 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
450 + NLKProductionCX + 443 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX > rChannel) {
451 isElastic = false; 444 isElastic = false;
452 isStrangeProduction = true; 445 isStrangeProduction = true;
453 // NN -> NNKKb channel is chosen 446 // NN -> NNKKb channel is chosen
454 INCL_DEBUG("NN interaction: NNKKb chan 447 INCL_DEBUG("NN interaction: NNKKb channel chosen" << '\n');
455 weight = bias_apply; 448 weight = bias_apply;
456 return new NNToNNKKbChannel(particle1, 449 return new NNToNNKKbChannel(particle1, particle2);
457 } else if(elasticCX + deltaProductionCX 450 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + fourPiProductionCX
458 + etaProductionCX + 451 + etaProductionCX + etadeltaProductionCX + etaonePiProductionCX + etatwoPiProductionCX + etathreePiProductionCX + etafourPiProductionCX
459 + omegaProductionCX 452 + omegaProductionCX + omegadeltaProductionCX + omegaonePiProductionCX + omegatwoPiProductionCX + omegathreePiProductionCX + omegafourPiProductionCX
460 + NLKProductionCX + 453 + NLKProductionCX + NLKpiProductionCX + NLK2piProductionCX + NSKProductionCX + NSKpiProductionCX + NSK2piProductionCX + NNKKbProductionCX + NNMissingCX> rChannel) {
461 isElastic = false; 454 isElastic = false;
462 isStrangeProduction = true; 455 isStrangeProduction = true;
463 // NN -> Missing Strangeness channel is chosen 456 // NN -> Missing Strangeness channel is chosen
464 INCL_DEBUG("NN interaction: Missing St 457 INCL_DEBUG("NN interaction: Missing Strangeness channel chosen" << '\n');
465 weight = bias_apply; 458 weight = bias_apply;
466 return new NNToMissingStrangenessChann 459 return new NNToMissingStrangenessChannel(particle1, particle2);
467 } else { 460 } else {
468 INCL_WARN("inconsistency within the NN 461 INCL_WARN("inconsistency within the NN Cross Sections (sum!=inelastic)" << '\n');
469 if(NNMissingCX>0.) { 462 if(NNMissingCX>0.) {
470 INCL_WARN("Returning an Missing St 463 INCL_WARN("Returning an Missing Strangeness channel" << '\n');
471 weight = bias_apply; 464 weight = bias_apply;
472 isElastic = false; 465 isElastic = false;
473 isStrangeProduction = true; 466 isStrangeProduction = true;
474 return new NNToNNKKbChannel(partic 467 return new NNToNNKKbChannel(particle1, particle2);
475 } else if(NNKKbProductionCX>0.) { 468 } else if(NNKKbProductionCX>0.) {
476 INCL_WARN("Returning an NNKKb chan 469 INCL_WARN("Returning an NNKKb channel" << '\n');
477 weight = bias_apply; 470 weight = bias_apply;
478 isElastic = false; 471 isElastic = false;
479 isStrangeProduction = true; 472 isStrangeProduction = true;
480 return new NNToNNKKbChannel(partic 473 return new NNToNNKKbChannel(particle1, particle2);
481 } else if(NSK2piProductionCX>0.) { 474 } else if(NSK2piProductionCX>0.) {
482 INCL_WARN("Returning an NSK2pi cha 475 INCL_WARN("Returning an NSK2pi channel" << '\n');
483 weight = bias_apply; 476 weight = bias_apply;
484 isElastic = false; 477 isElastic = false;
485 isStrangeProduction = true; 478 isStrangeProduction = true;
486 return new NNToNSK2piChannel(parti 479 return new NNToNSK2piChannel(particle1, particle2);
487 } else if(NSKpiProductionCX>0.) { 480 } else if(NSKpiProductionCX>0.) {
488 INCL_WARN("Returning an NSKpi chan 481 INCL_WARN("Returning an NSKpi channel" << '\n');
489 weight = bias_apply; 482 weight = bias_apply;
490 isElastic = false; 483 isElastic = false;
491 isStrangeProduction = true; 484 isStrangeProduction = true;
492 return new NNToNSKpiChannel(partic 485 return new NNToNSKpiChannel(particle1, particle2);
493 } else if(NSKProductionCX>0.) { 486 } else if(NSKProductionCX>0.) {
494 INCL_WARN("Returning an NSK channe 487 INCL_WARN("Returning an NSK channel" << '\n');
495 weight = bias_apply; 488 weight = bias_apply;
496 isElastic = false; 489 isElastic = false;
497 isStrangeProduction = true; 490 isStrangeProduction = true;
498 return new NNToNSKChannel(particle 491 return new NNToNSKChannel(particle1, particle2);
499 } else if(NLK2piProductionCX>0.) { 492 } else if(NLK2piProductionCX>0.) {
500 INCL_WARN("Returning an NLK2pi cha 493 INCL_WARN("Returning an NLK2pi channel" << '\n');
501 weight = bias_apply; 494 weight = bias_apply;
502 isElastic = false; 495 isElastic = false;
503 isStrangeProduction = true; 496 isStrangeProduction = true;
504 return new NNToNLK2piChannel(parti 497 return new NNToNLK2piChannel(particle1, particle2);
505 } else if(NLKpiProductionCX>0.) { 498 } else if(NLKpiProductionCX>0.) {
506 INCL_WARN("Returning an NLKpi chan 499 INCL_WARN("Returning an NLKpi channel" << '\n');
507 weight = bias_apply; 500 weight = bias_apply;
508 isElastic = false; 501 isElastic = false;
509 isStrangeProduction = true; 502 isStrangeProduction = true;
510 return new NNToNLKpiChannel(partic 503 return new NNToNLKpiChannel(particle1, particle2);
511 } else if(NLKProductionCX>0.) { 504 } else if(NLKProductionCX>0.) {
512 INCL_WARN("Returning an NLK channe 505 INCL_WARN("Returning an NLK channel" << '\n');
513 weight = bias_apply; 506 weight = bias_apply;
514 isElastic = false; 507 isElastic = false;
515 isStrangeProduction = true; 508 isStrangeProduction = true;
516 return new NNToNLKChannel(particle 509 return new NNToNLKChannel(particle1, particle2);
517 } else if(omegafourPiProductionCX>0.) 510 } else if(omegafourPiProductionCX>0.) {
518 INCL_WARN("Returning an Omega + fo 511 INCL_WARN("Returning an Omega + four Pions channel" << '\n');
519 weight = counterweight; 512 weight = counterweight;
520 isElastic = false; 513 isElastic = false;
521 return new NNOmegaToMultiPionsChan 514 return new NNOmegaToMultiPionsChannel(4,particle1, particle2);
522 } else if(omegathreePiProductionCX>0.) 515 } else if(omegathreePiProductionCX>0.) {
523 INCL_WARN("Returning an Omega + th 516 INCL_WARN("Returning an Omega + three Pions channel" << '\n');
524 weight = counterweight; 517 weight = counterweight;
525 isElastic = false; 518 isElastic = false;
526 return new NNOmegaToMultiPionsChan 519 return new NNOmegaToMultiPionsChannel(3,particle1, particle2);
527 } else if(omegatwoPiProductionCX>0.) { 520 } else if(omegatwoPiProductionCX>0.) {
528 INCL_WARN("Returning an Omega + tw 521 INCL_WARN("Returning an Omega + two Pions channel" << '\n');
529 weight = counterweight; 522 weight = counterweight;
530 isElastic = false; 523 isElastic = false;
531 return new NNOmegaToMultiPionsChan 524 return new NNOmegaToMultiPionsChannel(2,particle1, particle2);
532 } else if(omegaonePiProductionCX>0.) { 525 } else if(omegaonePiProductionCX>0.) {
533 INCL_WARN("Returning an Omega + on 526 INCL_WARN("Returning an Omega + one Pion channel" << '\n');
534 weight = counterweight; 527 weight = counterweight;
535 isElastic = false; 528 isElastic = false;
536 return new NNOmegaToMultiPionsChannel 529 return new NNOmegaToMultiPionsChannel(1,particle1, particle2);
537 } else if(omegadeltaProductionCX>0.) { 530 } else if(omegadeltaProductionCX>0.) {
538 INCL_WARN("Returning an Omega + De 531 INCL_WARN("Returning an Omega + Delta channel" << '\n');
539 weight = counterweight; 532 weight = counterweight;
540 isElastic = false; 533 isElastic = false;
541 return new NDeltaOmegaProductionCh 534 return new NDeltaOmegaProductionChannel(particle1, particle2);
542 } else if(omegaProductionCX>0.) { 535 } else if(omegaProductionCX>0.) {
543 INCL_WARN("Returning an Omega chan 536 INCL_WARN("Returning an Omega channel" << '\n');
544 weight = counterweight; 537 weight = counterweight;
545 isElastic = false; 538 isElastic = false;
546 return new NNToNNOmegaChannel(part 539 return new NNToNNOmegaChannel(particle1, particle2);
547 } else if(etafourPiProductionCX>0.) { 540 } else if(etafourPiProductionCX>0.) {
548 INCL_WARN("Returning an Eta + four 541 INCL_WARN("Returning an Eta + four Pions channel" << '\n');
549 weight = counterweight; 542 weight = counterweight;
550 isElastic = false; 543 isElastic = false;
551 return new NNEtaToMultiPionsChanne 544 return new NNEtaToMultiPionsChannel(4,particle1, particle2);
552 } else if(etathreePiProductionCX>0.) { 545 } else if(etathreePiProductionCX>0.) {
553 INCL_WARN("Returning an Eta + thre 546 INCL_WARN("Returning an Eta + threev channel" << '\n');
554 weight = counterweight; 547 weight = counterweight;
555 isElastic = false; 548 isElastic = false;
556 return new NNEtaToMultiPionsChanne 549 return new NNEtaToMultiPionsChannel(3,particle1, particle2);
557 } else if(etatwoPiProductionCX>0.) { 550 } else if(etatwoPiProductionCX>0.) {
558 INCL_WARN("Returning an Eta + two 551 INCL_WARN("Returning an Eta + two Pions channel" << '\n');
559 weight = counterweight; 552 weight = counterweight;
560 isElastic = false; 553 isElastic = false;
561 return new NNEtaToMultiPionsChanne 554 return new NNEtaToMultiPionsChannel(2,particle1, particle2);
562 } else if(etaonePiProductionCX>0.) { 555 } else if(etaonePiProductionCX>0.) {
563 INCL_WARN("Returning an Eta + one 556 INCL_WARN("Returning an Eta + one Pion channel" << '\n');
564 weight = counterweight; 557 weight = counterweight;
565 isElastic = false; 558 isElastic = false;
566 return new NNEtaToMultiPionsChanne 559 return new NNEtaToMultiPionsChannel(1,particle1, particle2);
567 } else if(etadeltaProductionCX>0.) { 560 } else if(etadeltaProductionCX>0.) {
568 INCL_WARN("Returning an Eta + Delt 561 INCL_WARN("Returning an Eta + Delta channel" << '\n');
569 weight = counterweight; 562 weight = counterweight;
570 isElastic = false; 563 isElastic = false;
571 return new NDeltaEtaProductionChan 564 return new NDeltaEtaProductionChannel(particle1, particle2);
572 } else if(etaProductionCX>0.) { 565 } else if(etaProductionCX>0.) {
573 INCL_WARN("Returning an Eta channe 566 INCL_WARN("Returning an Eta channel" << '\n');
574 weight = counterweight; 567 weight = counterweight;
575 isElastic = false; 568 isElastic = false;
576 return new NNToNNEtaChannel(partic 569 return new NNToNNEtaChannel(particle1, particle2);
577 } else if(fourPiProductionCX>0.) { 570 } else if(fourPiProductionCX>0.) {
578 INCL_WARN("Returning a 4pi channel 571 INCL_WARN("Returning a 4pi channel" << '\n');
579 weight = counterweight; 572 weight = counterweight;
580 isElastic = false; 573 isElastic = false;
581 return new NNToMultiPionsChannel(4 574 return new NNToMultiPionsChannel(4,particle1, particle2);
582 } else if(threePiProductionCX>0.) { 575 } else if(threePiProductionCX>0.) {
583 INCL_WARN("Returning a 3pi channel 576 INCL_WARN("Returning a 3pi channel" << '\n');
584 weight = counterweight; 577 weight = counterweight;
585 isElastic = false; 578 isElastic = false;
586 return new NNToMultiPionsChannel(3 579 return new NNToMultiPionsChannel(3,particle1, particle2);
587 } else if(twoPiProductionCX>0.) { 580 } else if(twoPiProductionCX>0.) {
588 INCL_WARN("Returning a 2pi channel 581 INCL_WARN("Returning a 2pi channel" << '\n');
589 weight = counterweight; 582 weight = counterweight;
590 isElastic = false; 583 isElastic = false;
591 return new NNToMultiPionsChannel(2 584 return new NNToMultiPionsChannel(2,particle1, particle2);
592 } else if(onePiProductionCX>0.) { 585 } else if(onePiProductionCX>0.) {
593 INCL_WARN("Returning a 1pi channel 586 INCL_WARN("Returning a 1pi channel" << '\n');
594 weight = counterweight; 587 weight = counterweight;
595 isElastic = false; 588 isElastic = false;
596 return new NNToMultiPionsChannel(1 589 return new NNToMultiPionsChannel(1,particle1, particle2);
597 } else if(deltaProductionCX>0.) { 590 } else if(deltaProductionCX>0.) {
598 INCL_WARN("Returning a delta-produ 591 INCL_WARN("Returning a delta-production channel" << '\n');
599 weight = counterweight; 592 weight = counterweight;
600 isElastic = false; 593 isElastic = false;
601 return new DeltaProductionChannel( 594 return new DeltaProductionChannel(particle1, particle2);
602 } else { 595 } else {
603 INCL_WARN("Returning an elastic ch 596 INCL_WARN("Returning an elastic channel" << '\n');
604 weight = counterweight; 597 weight = counterweight;
605 isElastic = true; 598 isElastic = true;
606 return new ElasticChannel(particle 599 return new ElasticChannel(particle1, particle2);
607 } 600 }
608 } 601 }
609 602
610 //// NDelta 603 //// NDelta
611 } 604 }
612 else if((particle1->isNucleon() && particl 605 else if((particle1->isNucleon() && particle2->isDelta()) ||
613 (particle1->isDelta() && part 606 (particle1->isDelta() && particle2->isNucleon())) {
614 607
615 G4double NLKProductionCX = CrossSect 608 G4double NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
616 G4double NSKProductionCX = CrossSect 609 G4double NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
617 G4double DeltaLKProductionCX = Cross 610 G4double DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
618 G4double DeltaSKProductionCX = Cross 611 G4double DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
619 G4double NNKKbProductionCX = CrossSe 612 G4double NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
620 613
621 const G4double UnStrangeProdCX = Cro 614 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::NDeltaToNN(particle1, particle2);
622 const G4double StrangenessProdCX = ( 615 const G4double StrangenessProdCX = (NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX)/bias_apply;
623 616
624 G4double counterweight = (1. - bias_ 617 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
625 618
626 if(counterweight < 0.5){ 619 if(counterweight < 0.5){
627 counterweight = 0.5; 620 counterweight = 0.5;
628 bias_apply = 0.5*UnStrangeProdCX/ 621 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
629 622
630 NLKProductionCX = CrossSections:: 623 NLKProductionCX = CrossSections::NDeltaToNLK(particle1, particle2)*bias_apply;
631 NSKProductionCX = CrossSections:: 624 NSKProductionCX = CrossSections::NDeltaToNSK(particle1, particle2)*bias_apply;
632 DeltaLKProductionCX = CrossSectio 625 DeltaLKProductionCX = CrossSections::NDeltaToDeltaLK(particle1, particle2)*bias_apply;
633 DeltaSKProductionCX = CrossSectio 626 DeltaSKProductionCX = CrossSections::NDeltaToDeltaSK(particle1, particle2)*bias_apply;
634 NNKKbProductionCX = CrossSections 627 NNKKbProductionCX = CrossSections::NDeltaToNNKKb(particle1, particle2)*bias_apply;
635 } 628 }
636 629
637 G4double elasticCX = CrossSections:: 630 G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
638 G4double recombinationCX = CrossSect 631 G4double recombinationCX = CrossSections::NDeltaToNN(particle1, particle2)*counterweight;
639 632
640 const G4double rChannel=Random::shoo 633 const G4double rChannel=Random::shoot() * (StrangenessProdCX + UnStrangeProdCX);
641 634
642 if(elasticCX > rChannel) { 635 if(elasticCX > rChannel) {
643 isElastic = true; 636 isElastic = true;
644 // Elastic N Delta channel 637 // Elastic N Delta channel
645 INCL_DEBUG("NDelta interaction: e 638 INCL_DEBUG("NDelta interaction: elastic channel chosen" << '\n');
646 weight = counterweight; 639 weight = counterweight;
647 return new ElasticChannel(particl 640 return new ElasticChannel(particle1, particle2);
648 } else if (elasticCX + recombination 641 } else if (elasticCX + recombinationCX > rChannel){
649 isElastic = false; 642 isElastic = false;
650 // Recombination 643 // Recombination
651 // NDelta -> NN channel is chosen 644 // NDelta -> NN channel is chosen
652 INCL_DEBUG("NDelta interaction: r 645 INCL_DEBUG("NDelta interaction: recombination channel chosen" << '\n');
653 weight = counterweight; 646 weight = counterweight;
654 return new RecombinationChannel(p 647 return new RecombinationChannel(particle1, particle2);
655 } else if (elasticCX + recombination 648 } else if (elasticCX + recombinationCX + NLKProductionCX > rChannel){
656 isElastic = false; 649 isElastic = false;
657 isStrangeProduction = true; 650 isStrangeProduction = true;
658 // NDelta -> NLK channel is chosen 651 // NDelta -> NLK channel is chosen
659 INCL_DEBUG("NDelta interaction: N 652 INCL_DEBUG("NDelta interaction: NLK channel chosen" << '\n');
660 weight = bias_apply; 653 weight = bias_apply;
661 return new NDeltaToNLKChannel(par 654 return new NDeltaToNLKChannel(particle1, particle2);
662 } else if (elasticCX + recombination 655 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX > rChannel){
663 isElastic = false; 656 isElastic = false;
664 isStrangeProduction = true; 657 isStrangeProduction = true;
665 // NDelta -> NSK channel is chosen 658 // NDelta -> NSK channel is chosen
666 INCL_DEBUG("NDelta interaction: N 659 INCL_DEBUG("NDelta interaction: NSK channel chosen" << '\n');
667 weight = bias_apply; 660 weight = bias_apply;
668 return new NDeltaToNSKChannel(par 661 return new NDeltaToNSKChannel(particle1, particle2);
669 } else if (elasticCX + recombination 662 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX > rChannel){
670 isElastic = false; 663 isElastic = false;
671 isStrangeProduction = true; 664 isStrangeProduction = true;
672 // NDelta -> DeltaLK channel is chosen 665 // NDelta -> DeltaLK channel is chosen
673 INCL_DEBUG("NDelta interaction: D 666 INCL_DEBUG("NDelta interaction: DeltaLK channel chosen" << '\n');
674 weight = bias_apply; 667 weight = bias_apply;
675 return new NDeltaToDeltaLKChannel 668 return new NDeltaToDeltaLKChannel(particle1, particle2);
676 } else if (elasticCX + recombination 669 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX > rChannel){
677 isElastic = false; 670 isElastic = false;
678 isStrangeProduction = true; 671 isStrangeProduction = true;
679 // NDelta -> DeltaSK channel is chosen 672 // NDelta -> DeltaSK channel is chosen
680 INCL_DEBUG("NDelta interaction: D 673 INCL_DEBUG("NDelta interaction: DeltaSK channel chosen" << '\n');
681 weight = bias_apply; 674 weight = bias_apply;
682 return new NDeltaToDeltaSKChannel 675 return new NDeltaToDeltaSKChannel(particle1, particle2);
683 } else if (elasticCX + recombination 676 } else if (elasticCX + recombinationCX + NLKProductionCX + NSKProductionCX + DeltaLKProductionCX + DeltaSKProductionCX + NNKKbProductionCX > rChannel){
684 isElastic = false; 677 isElastic = false;
685 isStrangeProduction = true; 678 isStrangeProduction = true;
686 // NDelta -> NNKKb channel is chosen 679 // NDelta -> NNKKb channel is chosen
687 INCL_DEBUG("NDelta interaction: N 680 INCL_DEBUG("NDelta interaction: NNKKb channel chosen" << '\n');
688 weight = bias_apply; 681 weight = bias_apply;
689 return new NDeltaToNNKKbChannel(p 682 return new NDeltaToNNKKbChannel(particle1, particle2);
690 } 683 }
691 else{ 684 else{
692 INCL_ERROR("rChannel > (Strangene 685 INCL_ERROR("rChannel > (StrangenessProdCX + UnStrangeProdCX) in NDelta interaction: return an elastic channel" << '\n');
693 weight = counterweight; 686 weight = counterweight;
694 isElastic = true; 687 isElastic = true;
695 return new ElasticChannel(particl 688 return new ElasticChannel(particle1, particle2);
696 } 689 }
697 690
698 //// DeltaDelta 691 //// DeltaDelta
699 } else if(particle1->isDelta() && particle 692 } else if(particle1->isDelta() && particle2->isDelta()) {
700 isElastic = true; 693 isElastic = true;
701 INCL_DEBUG("DeltaDelta interaction: el 694 INCL_DEBUG("DeltaDelta interaction: elastic channel chosen" << '\n');
702 return new ElasticChannel(particle1, p 695 return new ElasticChannel(particle1, particle2);
703 696
704 //// PiN 697 //// PiN
705 } else if(isPiN) { 698 } else if(isPiN) {
706 699
707 G4double LKProdCX = CrossSections::NpiTo 700 G4double LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
708 G4double SKProdCX = CrossSections::NpiTo 701 G4double SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
709 G4double LKpiProdCX = CrossSections::Npi 702 G4double LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
710 G4double SKpiProdCX = CrossSections::Npi 703 G4double SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
711 G4double LK2piProdCX = CrossSections::Np 704 G4double LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
712 G4double SK2piProdCX = CrossSections::Np 705 G4double SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
713 G4double NKKbProdCX = CrossSections::Npi 706 G4double NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
714 G4double MissingCX = CrossSections::NpiT 707 G4double MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
715 708
716 const G4double UnStrangeProdCX = CrossSe 709 const G4double UnStrangeProdCX = CrossSections::elastic(particle1, particle2) + CrossSections::piNToDelta(particle1, particle2)
717 + CrossSect 710 + CrossSections::piNToxPiN(2,particle1, particle2) + CrossSections::piNToxPiN(3,particle1, particle2) + CrossSections::piNToxPiN(4,particle1, particle2)
718 + CrossSect 711 + CrossSections::piNToEtaN(particle1, particle2) + CrossSections::piNToOmegaN(particle1, particle2);
719 const G4double StrangenessProdCX = (LKPr 712 const G4double StrangenessProdCX = (LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX)/bias_apply;
720 713
721 G4double counterweight = (1. - bias_appl 714 G4double counterweight = (1. - bias_apply * StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX))/(1. - StrangenessProdCX / (StrangenessProdCX + UnStrangeProdCX));
722 715
723 if(counterweight < 0.5) { 716 if(counterweight < 0.5) {
724 counterweight = 0.5; 717 counterweight = 0.5;
725 bias_apply = 0.5*UnStrangeProdCX/Stra 718 bias_apply = 0.5*UnStrangeProdCX/StrangenessProdCX+1;
726 LKProdCX = CrossSections::NpiToLK(par 719 LKProdCX = CrossSections::NpiToLK(particle1,particle2)*bias_apply;
727 SKProdCX = CrossSections::NpiToSK(par 720 SKProdCX = CrossSections::NpiToSK(particle1,particle2)*bias_apply;
728 LKpiProdCX = CrossSections::NpiToLKpi 721 LKpiProdCX = CrossSections::NpiToLKpi(particle1,particle2)*bias_apply;
729 SKpiProdCX = CrossSections::NpiToSKpi 722 SKpiProdCX = CrossSections::NpiToSKpi(particle1,particle2)*bias_apply;
730 LK2piProdCX = CrossSections::NpiToLK2 723 LK2piProdCX = CrossSections::NpiToLK2pi(particle1,particle2)*bias_apply;
731 SK2piProdCX = CrossSections::NpiToSK2 724 SK2piProdCX = CrossSections::NpiToSK2pi(particle1,particle2)*bias_apply;
732 NKKbProdCX = CrossSections::NpiToNKKb 725 NKKbProdCX = CrossSections::NpiToNKKb(particle1,particle2)*bias_apply;
733 MissingCX = CrossSections::NpiToMissi 726 MissingCX = CrossSections::NpiToMissingStrangeness(particle1,particle2)*bias_apply;
734 } 727 }
735 728
736 729
737 const G4double elasticCX = CrossSections 730 const G4double elasticCX = CrossSections::elastic(particle1, particle2)*counterweight;
738 const G4double deltaProductionCX = Cross 731 const G4double deltaProductionCX = CrossSections::piNToDelta(particle1, particle2)*counterweight;
739 const G4double onePiProductionCX = Cross 732 const G4double onePiProductionCX = CrossSections::piNToxPiN(2,particle1, particle2)*counterweight;
740 const G4double twoPiProductionCX = Cross 733 const G4double twoPiProductionCX = CrossSections::piNToxPiN(3,particle1, particle2)*counterweight;
741 const G4double threePiProductionCX = 734 const G4double threePiProductionCX = CrossSections::piNToxPiN(4,particle1, particle2)*counterweight;
742 const G4double etaProductionCX = Cros 735 const G4double etaProductionCX = CrossSections::piNToEtaN(particle1, particle2)*counterweight;
743 const G4double omegaProductionCX = Cr 736 const G4double omegaProductionCX = CrossSections::piNToOmegaN(particle1, particle2)*counterweight;
744 737
745 const G4double totCX=CrossSections::tota 738 const G4double totCX=CrossSections::total(particle1, particle2);
746 739
747 // assert(std::fabs(totCX-elasticCX-deltaProdu 740 // assert(std::fabs(totCX-elasticCX-deltaProductionCX-onePiProductionCX-twoPiProductionCX-threePiProductionCX-etaProductionCX-omegaProductionCX-LKProdCX-SKProdCX-LKpiProdCX-SKpiProdCX-LK2piProdCX-SK2piProdCX-NKKbProdCX-MissingCX) < 0.15);
748 741
749 const G4double rChannel=Random::shoot() 742 const G4double rChannel=Random::shoot() * totCX;
750 743
751 if(elasticCX > rChannel) { 744 if(elasticCX > rChannel) {
752 isElastic = true; 745 isElastic = true;
753 // Elastic PiN channel 746 // Elastic PiN channel
754 INCL_DEBUG("PiN interaction: elastic c 747 INCL_DEBUG("PiN interaction: elastic channel chosen" << '\n');
755 weight = counterweight; 748 weight = counterweight;
756 return new PiNElasticChannel(particle1 749 return new PiNElasticChannel(particle1, particle2);
757 } else if(elasticCX + deltaProductionCX 750 } else if(elasticCX + deltaProductionCX > rChannel) {
758 isElastic = false; 751 isElastic = false;
759 // PiN -> Delta channel is chosen 752 // PiN -> Delta channel is chosen
760 INCL_DEBUG("PiN interaction: Delta cha 753 INCL_DEBUG("PiN interaction: Delta channel chosen" << '\n');
761 weight = counterweight; 754 weight = counterweight;
762 return new PiNToDeltaChannel(particle1 755 return new PiNToDeltaChannel(particle1, particle2);
763 } else if(elasticCX + deltaProductionCX 756 } else if(elasticCX + deltaProductionCX + onePiProductionCX > rChannel) {
764 isElastic = false; 757 isElastic = false;
765 // PiN -> PiNPi channel is chosen 758 // PiN -> PiNPi channel is chosen
766 INCL_DEBUG("PiN interaction: one Pion 759 INCL_DEBUG("PiN interaction: one Pion channel chosen" << '\n');
767 weight = counterweight; 760 weight = counterweight;
768 return new PiNToMultiPionsChannel(2,pa 761 return new PiNToMultiPionsChannel(2,particle1, particle2);
769 } else if(elasticCX + deltaProductionCX 762 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX > rChannel) {
770 isElastic = false; 763 isElastic = false;
771 // PiN -> PiN2Pi channel is chosen 764 // PiN -> PiN2Pi channel is chosen
772 INCL_DEBUG("PiN interaction: two Pions 765 INCL_DEBUG("PiN interaction: two Pions channel chosen" << '\n');
773 weight = counterweight; 766 weight = counterweight;
774 return new PiNToMultiPionsChannel(3,pa 767 return new PiNToMultiPionsChannel(3,particle1, particle2);
775 } else if(elasticCX + deltaProductionCX 768 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX > rChannel) {
776 isElastic = false; 769 isElastic = false;
777 // PiN -> PiN3Pi channel is chosen 770 // PiN -> PiN3Pi channel is chosen
778 INCL_DEBUG("PiN interaction: three Pio 771 INCL_DEBUG("PiN interaction: three Pions channel chosen" << '\n');
779 weight = counterweight; 772 weight = counterweight;
780 return new PiNToMultiPionsChannel(4,pa 773 return new PiNToMultiPionsChannel(4,particle1, particle2);
781 } else if(elasticCX + deltaProductionCX 774 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX > rChannel) {
782 isElastic = false; 775 isElastic = false;
783 // PiN -> EtaN channel is chosen 776 // PiN -> EtaN channel is chosen
784 INCL_DEBUG("PiN interaction: Eta chann 777 INCL_DEBUG("PiN interaction: Eta channel chosen" << '\n');
785 weight = counterweight; 778 weight = counterweight;
786 return new PiNToEtaChannel(particle1, 779 return new PiNToEtaChannel(particle1, particle2);
787 } else if(elasticCX + deltaProductionCX 780 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX > rChannel) {
788 isElastic = false; 781 isElastic = false;
789 // PiN -> OmegaN channel is chosen 782 // PiN -> OmegaN channel is chosen
790 INCL_DEBUG("PiN interaction: Omega cha 783 INCL_DEBUG("PiN interaction: Omega channel chosen" << '\n');
791 weight = counterweight; 784 weight = counterweight;
792 return new PiNToOmegaChannel(particle1 785 return new PiNToOmegaChannel(particle1, particle2);
793 } else if(elasticCX + deltaProductionCX 786 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
794 + LKProdCX > rChanne 787 + LKProdCX > rChannel) {
795 isElastic = false; 788 isElastic = false;
796 isStrangeProduction = true; 789 isStrangeProduction = true;
797 // PiN -> LK channel is chosen 790 // PiN -> LK channel is chosen
798 INCL_DEBUG("PiN interaction: LK channe 791 INCL_DEBUG("PiN interaction: LK channel chosen" << '\n');
799 weight = bias_apply; 792 weight = bias_apply;
800 return new NpiToLKChannel(particle1, p 793 return new NpiToLKChannel(particle1, particle2);
801 } else if(elasticCX + deltaProductionCX 794 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
802 + LKProdCX + SKProdC 795 + LKProdCX + SKProdCX > rChannel) {
803 isElastic = false; 796 isElastic = false;
804 isStrangeProduction = true; 797 isStrangeProduction = true;
805 // PiN -> SK channel is chosen 798 // PiN -> SK channel is chosen
806 INCL_DEBUG("PiN interaction: SK channe 799 INCL_DEBUG("PiN interaction: SK channel chosen" << '\n');
807 weight = bias_apply; 800 weight = bias_apply;
808 return new NpiToSKChannel(particle1, p 801 return new NpiToSKChannel(particle1, particle2);
809 } else if(elasticCX + deltaProductionCX 802 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
810 + LKProdCX + SKProdC 803 + LKProdCX + SKProdCX + LKpiProdCX > rChannel) {
811 isElastic = false; 804 isElastic = false;
812 isStrangeProduction = true; 805 isStrangeProduction = true;
813 // PiN -> LKpi channel is chosen 806 // PiN -> LKpi channel is chosen
814 INCL_DEBUG("PiN interaction: LKpi chan 807 INCL_DEBUG("PiN interaction: LKpi channel chosen" << '\n');
815 weight = bias_apply; 808 weight = bias_apply;
816 return new NpiToLKpiChannel(particle1, 809 return new NpiToLKpiChannel(particle1, particle2);
817 } else if(elasticCX + deltaProductionCX 810 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
818 + LKProdCX + SKProdC 811 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX > rChannel) {
819 isElastic = false; 812 isElastic = false;
820 isStrangeProduction = true; 813 isStrangeProduction = true;
821 // PiN -> SKpi channel is chosen 814 // PiN -> SKpi channel is chosen
822 INCL_DEBUG("PiN interaction: SKpi chan 815 INCL_DEBUG("PiN interaction: SKpi channel chosen" << '\n');
823 weight = bias_apply; 816 weight = bias_apply;
824 return new NpiToSKpiChannel(particle1, 817 return new NpiToSKpiChannel(particle1, particle2);
825 } else if(elasticCX + deltaProductionCX 818 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
826 + LKProdCX + SKProdC 819 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX > rChannel) {
827 isElastic = false; 820 isElastic = false;
828 isStrangeProduction = true; 821 isStrangeProduction = true;
829 // PiN -> LK2pi channel is chosen 822 // PiN -> LK2pi channel is chosen
830 INCL_DEBUG("PiN interaction: LK2pi cha 823 INCL_DEBUG("PiN interaction: LK2pi channel chosen" << '\n');
831 weight = bias_apply; 824 weight = bias_apply;
832 return new NpiToLK2piChannel(particle1 825 return new NpiToLK2piChannel(particle1, particle2);
833 } else if(elasticCX + deltaProductionCX 826 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
834 + LKProdCX + SKProdC 827 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX > rChannel) {
835 isElastic = false; 828 isElastic = false;
836 isStrangeProduction = true; 829 isStrangeProduction = true;
837 // PiN -> SK2pi channel is chosen 830 // PiN -> SK2pi channel is chosen
838 INCL_DEBUG("PiN interaction: SK2pi cha 831 INCL_DEBUG("PiN interaction: SK2pi channel chosen" << '\n');
839 weight = bias_apply; 832 weight = bias_apply;
840 return new NpiToSK2piChannel(particle1 833 return new NpiToSK2piChannel(particle1, particle2);
841 } else if(elasticCX + deltaProductionCX 834 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
842 + LKProdCX + SKProdC 835 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX > rChannel) {
843 isElastic = false; 836 isElastic = false;
844 isStrangeProduction = true; 837 isStrangeProduction = true;
845 // PiN -> NKKb channel is chosen 838 // PiN -> NKKb channel is chosen
846 INCL_DEBUG("PiN interaction: NKKb chan 839 INCL_DEBUG("PiN interaction: NKKb channel chosen" << '\n');
847 weight = bias_apply; 840 weight = bias_apply;
848 return new NpiToNKKbChannel(particle1, 841 return new NpiToNKKbChannel(particle1, particle2);
849 } else if(elasticCX + deltaProductionCX 842 } else if(elasticCX + deltaProductionCX + onePiProductionCX + twoPiProductionCX + threePiProductionCX + etaProductionCX+ omegaProductionCX
850 + LKProdCX + SKProdC 843 + LKProdCX + SKProdCX + LKpiProdCX + SKpiProdCX + LK2piProdCX + SK2piProdCX + NKKbProdCX + MissingCX> rChannel) {
851 isElastic = false; 844 isElastic = false;
852 isStrangeProduction = true; 845 isStrangeProduction = true;
853 // PiN -> Missinge Strangeness channel is chos 846 // PiN -> Missinge Strangeness channel is chosen
854 INCL_DEBUG("PiN interaction: Missinge 847 INCL_DEBUG("PiN interaction: Missinge Strangeness channel chosen" << '\n');
855 weight = bias_apply; 848 weight = bias_apply;
856 return new NpiToMissingStrangenessChan 849 return new NpiToMissingStrangenessChannel(particle1, particle2);
857 } 850 }
858 else { 851 else {
859 INCL_WARN("inconsistency within the P 852 INCL_WARN("inconsistency within the PiN Cross Sections (sum!=inelastic)" << '\n');
860 if(MissingCX>0.) { 853 if(MissingCX>0.) {
861 INCL_WARN("Returning a Missinge St 854 INCL_WARN("Returning a Missinge Strangeness channel" << '\n');
862 weight = bias_apply; 855 weight = bias_apply;
863 isElastic = false; 856 isElastic = false;
864 isStrangeProduction = true; 857 isStrangeProduction = true;
865 return new NpiToMissingStrangeness 858 return new NpiToMissingStrangenessChannel(particle1, particle2);
866 } else if(NKKbProdCX>0.) { 859 } else if(NKKbProdCX>0.) {
867 INCL_WARN("Returning a NKKb channe 860 INCL_WARN("Returning a NKKb channel" << '\n');
868 weight = bias_apply; 861 weight = bias_apply;
869 isElastic = false; 862 isElastic = false;
870 isStrangeProduction = true; 863 isStrangeProduction = true;
871 return new NpiToNKKbChannel(partic 864 return new NpiToNKKbChannel(particle1, particle2);
872 } else if(SK2piProdCX>0.) { 865 } else if(SK2piProdCX>0.) {
873 INCL_WARN("Returning a SK2pi chann 866 INCL_WARN("Returning a SK2pi channel" << '\n');
874 weight = bias_apply; 867 weight = bias_apply;
875 isElastic = false; 868 isElastic = false;
876 isStrangeProduction = true; 869 isStrangeProduction = true;
877 return new NpiToSK2piChannel(parti 870 return new NpiToSK2piChannel(particle1, particle2);
878 } else if(LK2piProdCX>0.) { 871 } else if(LK2piProdCX>0.) {
879 INCL_WARN("Returning a LK2pi chann 872 INCL_WARN("Returning a LK2pi channel" << '\n');
880 weight = bias_apply; 873 weight = bias_apply;
881 isElastic = false; 874 isElastic = false;
882 isStrangeProduction = true; 875 isStrangeProduction = true;
883 return new NpiToLK2piChannel(parti 876 return new NpiToLK2piChannel(particle1, particle2);
884 } else if(SKpiProdCX>0.) { 877 } else if(SKpiProdCX>0.) {
885 INCL_WARN("Returning a SKpi channe 878 INCL_WARN("Returning a SKpi channel" << '\n');
886 weight = bias_apply; 879 weight = bias_apply;
887 isElastic = false; 880 isElastic = false;
888 isStrangeProduction = true; 881 isStrangeProduction = true;
889 return new NpiToSKpiChannel(partic 882 return new NpiToSKpiChannel(particle1, particle2);
890 } else if(LKpiProdCX>0.) { 883 } else if(LKpiProdCX>0.) {
891 INCL_WARN("Returning a LKpi channe 884 INCL_WARN("Returning a LKpi channel" << '\n');
892 weight = bias_apply; 885 weight = bias_apply;
893 isElastic = false; 886 isElastic = false;
894 isStrangeProduction = true; 887 isStrangeProduction = true;
895 return new NpiToLKpiChannel(partic 888 return new NpiToLKpiChannel(particle1, particle2);
896 } else if(SKProdCX>0.) { 889 } else if(SKProdCX>0.) {
897 INCL_WARN("Returning a SK channel" 890 INCL_WARN("Returning a SK channel" << '\n');
898 weight = bias_apply; 891 weight = bias_apply;
899 isElastic = false; 892 isElastic = false;
900 isStrangeProduction = true; 893 isStrangeProduction = true;
901 return new NpiToSKChannel(particle 894 return new NpiToSKChannel(particle1, particle2);
902 } else if(LKProdCX>0.) { 895 } else if(LKProdCX>0.) {
903 INCL_WARN("Returning a LK channel" 896 INCL_WARN("Returning a LK channel" << '\n');
904 weight = bias_apply; 897 weight = bias_apply;
905 isElastic = false; 898 isElastic = false;
906 isStrangeProduction = true; 899 isStrangeProduction = true;
907 return new NpiToLKChannel(particle 900 return new NpiToLKChannel(particle1, particle2);
908 } else if(omegaProductionCX>0.) { 901 } else if(omegaProductionCX>0.) {
909 INCL_WARN("Returning a Omega chann 902 INCL_WARN("Returning a Omega channel" << '\n');
910 weight = counterweight; 903 weight = counterweight;
911 isElastic = false; 904 isElastic = false;
912 return new PiNToOmegaChannel(parti 905 return new PiNToOmegaChannel(particle1, particle2);
913 } else if(etaProductionCX>0.) { 906 } else if(etaProductionCX>0.) {
914 INCL_WARN("Returning a Eta channel 907 INCL_WARN("Returning a Eta channel" << '\n');
915 weight = counterweight; 908 weight = counterweight;
916 isElastic = false; 909 isElastic = false;
917 return new PiNToEtaChannel(particl 910 return new PiNToEtaChannel(particle1, particle2);
918 } else if(threePiProductionCX>0.) { 911 } else if(threePiProductionCX>0.) {
919 INCL_WARN("Returning a 3pi channel 912 INCL_WARN("Returning a 3pi channel" << '\n');
920 weight = counterweight; 913 weight = counterweight;
921 isElastic = false; 914 isElastic = false;
922 return new PiNToMultiPionsChannel( 915 return new PiNToMultiPionsChannel(4,particle1, particle2);
923 } else if(twoPiProductionCX>0.) { 916 } else if(twoPiProductionCX>0.) {
924 INCL_WARN("Returning a 2pi channel 917 INCL_WARN("Returning a 2pi channel" << '\n');
925 weight = counterweight; 918 weight = counterweight;
926 isElastic = false; 919 isElastic = false;
927 return new PiNToMultiPionsChannel( 920 return new PiNToMultiPionsChannel(3,particle1, particle2);
928 } else if(onePiProductionCX>0.) { 921 } else if(onePiProductionCX>0.) {
929 INCL_WARN("Returning a 1pi channel 922 INCL_WARN("Returning a 1pi channel" << '\n');
930 weight = counterweight; 923 weight = counterweight;
931 isElastic = false; 924 isElastic = false;
932 return new PiNToMultiPionsChannel( 925 return new PiNToMultiPionsChannel(2,particle1, particle2);
933 } else if(deltaProductionCX>0.) { 926 } else if(deltaProductionCX>0.) {
934 INCL_WARN("Returning a delta-produ 927 INCL_WARN("Returning a delta-production channel" << '\n');
935 weight = counterweight; 928 weight = counterweight;
936 isElastic = false; 929 isElastic = false;
937 return new PiNToDeltaChannel(parti 930 return new PiNToDeltaChannel(particle1, particle2);
938 } else { 931 } else {
939 INCL_WARN("Returning an elastic ch 932 INCL_WARN("Returning an elastic channel" << '\n');
940 weight = counterweight; 933 weight = counterweight;
941 isElastic = true; 934 isElastic = true;
942 return new PiNElasticChannel(parti 935 return new PiNElasticChannel(particle1, particle2);
943 } 936 }
944 } 937 }
945 } else if ((particle1->isNucleon() && part 938 } else if ((particle1->isNucleon() && particle2->isEta()) || (particle2->isNucleon() && particle1->isEta())) {
946 //// EtaN 939 //// EtaN
947 940
948 const G4double elasticCX = 941 const G4double elasticCX = CrossSections::elastic(particle1, particle2);
949 const G4double onePiProduc 942 const G4double onePiProductionCX = CrossSections::etaNToPiN(particle1, particle2);
950 const G4double twoPiProduc 943 const G4double twoPiProductionCX = CrossSections::etaNToPiPiN(particle1, particle2);
951 const G4double totCX=Cross 944 const G4double totCX=CrossSections::total(particle1, particle2);
952 // assert(std::fabs(totCX-elasticCX-onePiProdu 945 // assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
953 946
954 const G4double rChannel=Ra 947 const G4double rChannel=Random::shoot() * totCX;
955 948
956 if(elasticCX > rChannel) { 949 if(elasticCX > rChannel) {
957 // Elastic EtaN channel 950 // Elastic EtaN channel
958 isElastic = true; 951 isElastic = true;
959 INCL_DEBUG("EtaN inter 952 INCL_DEBUG("EtaN interaction: elastic channel chosen" << '\n');
960 return new EtaNElastic 953 return new EtaNElasticChannel(particle1, particle2);
961 } else if(elasticCX + oneP 954 } else if(elasticCX + onePiProductionCX > rChannel) {
962 isElastic = false; 955 isElastic = false;
963 // EtaN -> EtaPiN channel is chosen 956 // EtaN -> EtaPiN channel is chosen
964 INCL_DEBUG("EtaN inter 957 INCL_DEBUG("EtaN interaction: PiN channel chosen" << '\n');
965 return new EtaNToPiNCh 958 return new EtaNToPiNChannel(particle1, particle2);
966 } else if(elasticCX + oneP 959 } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
967 isElastic = false; 960 isElastic = false;
968 // EtaN -> EtaPiPiN channel is chosen 961 // EtaN -> EtaPiPiN channel is chosen
969 INCL_DEBUG("EtaN inter 962 INCL_DEBUG("EtaN interaction: PiPiN channel chosen" << '\n');
970 return new EtaNToPiPiN 963 return new EtaNToPiPiNChannel(particle1, particle2);
971 } 964 }
972 965
973 else { 966 else {
974 INCL_WARN("inconsisten 967 INCL_WARN("inconsistency within the EtaN Cross Sections (sum!=inelastic)" << '\n');
975 if(twoPiProductionCX>0 968 if(twoPiProductionCX>0.) {
976 INCL_WARN("Returni 969 INCL_WARN("Returning a PiPiN channel" << '\n');
977 isElastic = false; 970 isElastic = false;
978 return new EtaNToP 971 return new EtaNToPiPiNChannel(particle1, particle2);
979 } else if(onePiProduct 972 } else if(onePiProductionCX>0.) {
980 INCL_WARN("Returni 973 INCL_WARN("Returning a PiN channel" << '\n');
981 isElastic = false; 974 isElastic = false;
982 return new EtaNToP 975 return new EtaNToPiNChannel(particle1, particle2);
983 } else { 976 } else {
984 INCL_WARN("Returni 977 INCL_WARN("Returning an elastic channel" << '\n');
985 isElastic = true; 978 isElastic = true;
986 return new EtaNEla 979 return new EtaNElasticChannel(particle1, particle2);
987 } 980 }
988 } 981 }
989 982
990 } else if ((particle1->isNucleon() && part 983 } else if ((particle1->isNucleon() && particle2->isOmega()) || (particle2->isNucleon() && particle1->isOmega())) {
991 //// OmegaN 984 //// OmegaN
992 985
993 const G4double elasticCX = 986 const G4double elasticCX = CrossSections::elastic(particle1, particle2);
994 const G4double onePiProduc 987 const G4double onePiProductionCX = CrossSections::omegaNToPiN(particle1, particle2);
995 const G4double twoPiProduc 988 const G4double twoPiProductionCX = CrossSections::omegaNToPiPiN(particle1, particle2);
996 const G4double totCX=Cross 989 const G4double totCX=CrossSections::total(particle1, particle2);
997 // assert(std::fabs(totCX-elasticCX-onePiProdu 990 // assert(std::fabs(totCX-elasticCX-onePiProductionCX-twoPiProductionCX)<1.);
998 991
999 const G4double rChannel=Ra 992 const G4double rChannel=Random::shoot() * totCX;
1000 993
1001 if(elasticCX > rChannel) 994 if(elasticCX > rChannel) {
1002 // Elastic OmegaN channel 995 // Elastic OmegaN channel
1003 isElastic = true; 996 isElastic = true;
1004 INCL_DEBUG("OmegaN in 997 INCL_DEBUG("OmegaN interaction: elastic channel chosen" << '\n');
1005 return new OmegaNElas 998 return new OmegaNElasticChannel(particle1, particle2);
1006 } else if(elasticCX + one 999 } else if(elasticCX + onePiProductionCX > rChannel) {
1007 isElastic = false; 1000 isElastic = false;
1008 // OmegaN -> PiN channel is chosen 1001 // OmegaN -> PiN channel is chosen
1009 INCL_DEBUG("OmegaN interaction: P 1002 INCL_DEBUG("OmegaN interaction: PiN channel chosen" << '\n');
1010 return new OmegaNToPiNChannel(par 1003 return new OmegaNToPiNChannel(particle1, particle2);
1011 } else if(elasticCX + one 1004 } else if(elasticCX + onePiProductionCX + twoPiProductionCX > rChannel) {
1012 isElastic = false; 1005 isElastic = false;
1013 // OmegaN -> PiPiN channel is chosen 1006 // OmegaN -> PiPiN channel is chosen
1014 INCL_DEBUG("OmegaN in 1007 INCL_DEBUG("OmegaN interaction: PiPiN channel chosen" << '\n');
1015 return new OmegaNToPi 1008 return new OmegaNToPiPiNChannel(particle1, particle2);
1016 } 1009 }
1017 else { 1010 else {
1018 INCL_WARN("inconsiste 1011 INCL_WARN("inconsistency within the OmegaN Cross Sections (sum!=inelastic)" << '\n');
1019 if(twoPiProductionCX> 1012 if(twoPiProductionCX>0.) {
1020 INCL_WARN("Return 1013 INCL_WARN("Returning a PiPiN channel" << '\n');
1021 isElastic = false 1014 isElastic = false;
1022 return new OmegaN 1015 return new OmegaNToPiPiNChannel(particle1, particle2);
1023 } else if(onePiProduc 1016 } else if(onePiProductionCX>0.) {
1024 INCL_WARN("Return 1017 INCL_WARN("Returning a PiN channel" << '\n');
1025 isElastic = false 1018 isElastic = false;
1026 return new OmegaN 1019 return new OmegaNToPiNChannel(particle1, particle2);
1027 } else { 1020 } else {
1028 INCL_WARN("Return 1021 INCL_WARN("Returning an elastic channel" << '\n');
1029 isElastic = true; 1022 isElastic = true;
1030 return new OmegaN 1023 return new OmegaNElasticChannel(particle1, particle2);
1031 } 1024 }
1032 } 1025 }
1033 } else if ((particle1->isNucleon() && par 1026 } else if ((particle1->isNucleon() && particle2->isKaon()) || (particle2->isNucleon() && particle1->isKaon())) {
1034 //// KN 1027 //// KN
1035 const G4double elasticCX = CrossSecti 1028 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1036 const G4double quasielasticCX = Cross 1029 const G4double quasielasticCX = CrossSections::NKToNK(particle1,particle2);
1037 const G4double NKToNKpiCX = CrossSect 1030 const G4double NKToNKpiCX = CrossSections::NKToNKpi(particle1,particle2);
1038 const G4double NKToNK2piCX = CrossSec 1031 const G4double NKToNK2piCX = CrossSections::NKToNK2pi(particle1,particle2);
1039 const G4double totCX=CrossSections::t 1032 const G4double totCX=CrossSections::total(particle1, particle2);
1040 // assert(std::fabs(totCX-elasticCX-quasielas 1033 // assert(std::fabs(totCX-elasticCX-quasielasticCX-NKToNKpiCX-NKToNK2piCX)<0.1);
1041 1034
1042 const G4double rChannel=Random::shoot 1035 const G4double rChannel=Random::shoot() * totCX;
1043 if(elasticCX > rChannel){ 1036 if(elasticCX > rChannel){
1044 // Elastic KN channel is chosen 1037 // Elastic KN channel is chosen
1045 isElastic = true; 1038 isElastic = true;
1046 INCL_DEBUG("KN interaction: elast 1039 INCL_DEBUG("KN interaction: elastic channel chosen" << '\n');
1047 return new NKElasticChannel(parti 1040 return new NKElasticChannel(particle1, particle2);
1048 } else if(elasticCX + quasielasticCX 1041 } else if(elasticCX + quasielasticCX > rChannel){
1049 // Quasi-elastic KN channel is chosen 1042 // Quasi-elastic KN channel is chosen
1050 isElastic = false; // true ?? 1043 isElastic = false; // true ??
1051 INCL_DEBUG("KN interaction: quasi 1044 INCL_DEBUG("KN interaction: quasi-elastic channel chosen" << '\n');
1052 return new NKToNKChannel(particle 1045 return new NKToNKChannel(particle1, particle2);
1053 } else if(elasticCX + quasielasticCX 1046 } else if(elasticCX + quasielasticCX + NKToNKpiCX > rChannel){
1054 // KN -> NKpi channel is chosen 1047 // KN -> NKpi channel is chosen
1055 isElastic = false; 1048 isElastic = false;
1056 INCL_DEBUG("KN interaction: NKpi 1049 INCL_DEBUG("KN interaction: NKpi channel chosen" << '\n');
1057 return new NKToNKpiChannel(partic 1050 return new NKToNKpiChannel(particle1, particle2);
1058 } else if(elasticCX + quasielasticCX 1051 } else if(elasticCX + quasielasticCX + NKToNKpiCX + NKToNK2piCX > rChannel){
1059 // KN -> NK2pi channel is chosen 1052 // KN -> NK2pi channel is chosen
1060 isElastic = false; 1053 isElastic = false;
1061 INCL_DEBUG("KN interaction: NK2pi 1054 INCL_DEBUG("KN interaction: NK2pi channel chosen" << '\n');
1062 return new NKToNK2piChannel(parti 1055 return new NKToNK2piChannel(particle1, particle2);
1063 } else { 1056 } else {
1064 INCL_WARN("inconsistency within t 1057 INCL_WARN("inconsistency within the KN Cross Sections (sum!=inelastic)" << '\n');
1065 if(NKToNK2piCX>0.) { 1058 if(NKToNK2piCX>0.) {
1066 INCL_WARN("Returning a NKToNK 1059 INCL_WARN("Returning a NKToNK2pi channel" << '\n');
1067 isElastic = false; 1060 isElastic = false;
1068 return new NKToNK2piChannel(p 1061 return new NKToNK2piChannel(particle1, particle2);
1069 } else if(NKToNKpiCX>0.) { 1062 } else if(NKToNKpiCX>0.) {
1070 INCL_WARN("Returning a NKToNK 1063 INCL_WARN("Returning a NKToNKpi channel" << '\n');
1071 isElastic = false; 1064 isElastic = false;
1072 return new NKToNKpiChannel(pa 1065 return new NKToNKpiChannel(particle1, particle2);
1073 } else if(quasielasticCX>0.) { 1066 } else if(quasielasticCX>0.) {
1074 INCL_WARN("Returning a quasi- 1067 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1075 isElastic = false; // true ?? 1068 isElastic = false; // true ??
1076 return new NKToNKChannel(part 1069 return new NKToNKChannel(particle1, particle2);
1077 } else { 1070 } else {
1078 INCL_WARN("Returning an elast 1071 INCL_WARN("Returning an elastic channel" << '\n');
1079 isElastic = true; 1072 isElastic = true;
1080 return new NKElasticChannel(p 1073 return new NKElasticChannel(particle1, particle2);
1081 } 1074 }
1082 } 1075 }
1083 } else if ((particle1->isNucleon() && par 1076 } else if ((particle1->isNucleon() && particle2->isAntiKaon()) || (particle2->isNucleon() && particle1->isAntiKaon())) {
1084 //// KbN 1077 //// KbN
1085 const G4double elasticCX = CrossSecti 1078 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1086 const G4double quasielasticCX = Cross 1079 const G4double quasielasticCX = CrossSections::NKbToNKb(particle1,particle2);
1087 const G4double NKbToNKbpiCX = CrossSe 1080 const G4double NKbToNKbpiCX = CrossSections::NKbToNKbpi(particle1,particle2);
1088 const G4double NKbToNKb2piCX = CrossS 1081 const G4double NKbToNKb2piCX = CrossSections::NKbToNKb2pi(particle1,particle2);
1089 const G4double NKbToLpiCX = CrossSect 1082 const G4double NKbToLpiCX = CrossSections::NKbToLpi(particle1,particle2);
1090 const G4double NKbToL2piCX = CrossSec 1083 const G4double NKbToL2piCX = CrossSections::NKbToL2pi(particle1,particle2);
1091 const G4double NKbToSpiCX = CrossSect 1084 const G4double NKbToSpiCX = CrossSections::NKbToSpi(particle1,particle2);
1092 const G4double NKbToS2piCX = CrossSec 1085 const G4double NKbToS2piCX = CrossSections::NKbToS2pi(particle1,particle2);
1093 const G4double totCX=CrossSections::t 1086 const G4double totCX=CrossSections::total(particle1, particle2);
1094 // assert(std::fabs(totCX-elasticCX-quasielas 1087 // assert(std::fabs(totCX-elasticCX-quasielasticCX-NKbToNKbpiCX-NKbToNKb2piCX-NKbToLpiCX-NKbToL2piCX-NKbToSpiCX-NKbToS2piCX)<0.1);
1095 1088
1096 const G4double rChannel=Random::shoot 1089 const G4double rChannel=Random::shoot() * totCX;
1097 if(elasticCX > rChannel){ 1090 if(elasticCX > rChannel){
1098 // Elastic KbN channel is chosen 1091 // Elastic KbN channel is chosen
1099 isElastic = true; 1092 isElastic = true;
1100 INCL_DEBUG("KbN interaction: elas 1093 INCL_DEBUG("KbN interaction: elastic channel chosen" << '\n');
1101 return new NKbElasticChannel(part 1094 return new NKbElasticChannel(particle1, particle2);
1102 } else if(elasticCX + quasielasticCX 1095 } else if(elasticCX + quasielasticCX > rChannel){
1103 // Quasi-elastic KbN channel is chosen 1096 // Quasi-elastic KbN channel is chosen
1104 isElastic = false; // true ?? 1097 isElastic = false; // true ??
1105 INCL_DEBUG("KbN interaction: quas 1098 INCL_DEBUG("KbN interaction: quasi-elastic channel chosen" << '\n');
1106 return new NKbToNKbChannel(partic 1099 return new NKbToNKbChannel(particle1, particle2);
1107 } else if(elasticCX + quasielasticCX 1100 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX > rChannel){
1108 // KbN -> NKbpi channel is chosen 1101 // KbN -> NKbpi channel is chosen
1109 isElastic = false; 1102 isElastic = false;
1110 INCL_DEBUG("KbN interaction: NKbp 1103 INCL_DEBUG("KbN interaction: NKbpi channel chosen" << '\n');
1111 return new NKbToNKbpiChannel(part 1104 return new NKbToNKbpiChannel(particle1, particle2);
1112 } else if(elasticCX + quasielasticCX 1105 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX > rChannel){
1113 // KbN -> NKb2pi channel is chosen 1106 // KbN -> NKb2pi channel is chosen
1114 isElastic = false; 1107 isElastic = false;
1115 INCL_DEBUG("KbN interaction: NKb2 1108 INCL_DEBUG("KbN interaction: NKb2pi channel chosen" << '\n');
1116 return new NKbToNKb2piChannel(par 1109 return new NKbToNKb2piChannel(particle1, particle2);
1117 } else if(elasticCX + quasielasticCX 1110 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX > rChannel){
1118 // KbN -> Lpi channel is chosen 1111 // KbN -> Lpi channel is chosen
1119 isElastic = false; 1112 isElastic = false;
1120 INCL_DEBUG("KbN interaction: Lpi 1113 INCL_DEBUG("KbN interaction: Lpi channel chosen" << '\n');
1121 return new NKbToLpiChannel(partic 1114 return new NKbToLpiChannel(particle1, particle2);
1122 } else if(elasticCX + quasielasticCX 1115 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX > rChannel){
1123 // KbN -> L2pi channel is chosen 1116 // KbN -> L2pi channel is chosen
1124 isElastic = false; 1117 isElastic = false;
1125 INCL_DEBUG("KbN interaction: L2pi 1118 INCL_DEBUG("KbN interaction: L2pi channel chosen" << '\n');
1126 return new NKbToL2piChannel(parti 1119 return new NKbToL2piChannel(particle1, particle2);
1127 } else if(elasticCX + quasielasticCX 1120 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX > rChannel){
1128 // KbN -> Spi channel is chosen 1121 // KbN -> Spi channel is chosen
1129 isElastic = false; 1122 isElastic = false;
1130 INCL_DEBUG("KbN interaction: Spi 1123 INCL_DEBUG("KbN interaction: Spi channel chosen" << '\n');
1131 return new NKbToSpiChannel(partic 1124 return new NKbToSpiChannel(particle1, particle2);
1132 } else if(elasticCX + quasielasticCX 1125 } else if(elasticCX + quasielasticCX + NKbToNKbpiCX + NKbToNKb2piCX + NKbToLpiCX + NKbToL2piCX + NKbToSpiCX + NKbToS2piCX > rChannel){
1133 // KbN -> S2pi channel is chosen 1126 // KbN -> S2pi channel is chosen
1134 isElastic = false; 1127 isElastic = false;
1135 INCL_DEBUG("KbN interaction: S2pi 1128 INCL_DEBUG("KbN interaction: S2pi channel chosen" << '\n');
1136 return new NKbToS2piChannel(parti 1129 return new NKbToS2piChannel(particle1, particle2);
1137 } else { 1130 } else {
1138 INCL_WARN("inconsistency within t 1131 INCL_WARN("inconsistency within the KbN Cross Sections (sum!=inelastic)" << '\n');
1139 if(NKbToS2piCX>0.) { 1132 if(NKbToS2piCX>0.) {
1140 INCL_WARN("Returning a NKbToS 1133 INCL_WARN("Returning a NKbToS2pi channel" << '\n');
1141 isElastic = false; 1134 isElastic = false;
1142 return new NKbToS2piChannel(p 1135 return new NKbToS2piChannel(particle1, particle2);
1143 } else if(NKbToSpiCX>0.) { 1136 } else if(NKbToSpiCX>0.) {
1144 INCL_WARN("Returning a NKbToS 1137 INCL_WARN("Returning a NKbToSpi channel" << '\n');
1145 isElastic = false; 1138 isElastic = false;
1146 return new NKbToSpiChannel(pa 1139 return new NKbToSpiChannel(particle1, particle2);
1147 } else if(NKbToL2piCX>0.) { 1140 } else if(NKbToL2piCX>0.) {
1148 INCL_WARN("Returning a NKbToL 1141 INCL_WARN("Returning a NKbToL2pi channel" << '\n');
1149 isElastic = false; 1142 isElastic = false;
1150 return new NKbToL2piChannel(p 1143 return new NKbToL2piChannel(particle1, particle2);
1151 } else if(NKbToLpiCX>0.) { 1144 } else if(NKbToLpiCX>0.) {
1152 INCL_WARN("Returning a NKbToL 1145 INCL_WARN("Returning a NKbToLpi channel" << '\n');
1153 isElastic = false; 1146 isElastic = false;
1154 return new NKbToLpiChannel(pa 1147 return new NKbToLpiChannel(particle1, particle2);
1155 } else if(NKbToNKb2piCX>0.) { 1148 } else if(NKbToNKb2piCX>0.) {
1156 INCL_WARN("Returning a NKbToN 1149 INCL_WARN("Returning a NKbToNKb2pi channel" << '\n');
1157 isElastic = false; 1150 isElastic = false;
1158 return new NKbToNKb2piChannel 1151 return new NKbToNKb2piChannel(particle1, particle2);
1159 } else if(NKbToNKbpiCX>0.) { 1152 } else if(NKbToNKbpiCX>0.) {
1160 INCL_WARN("Returning a NKbToN 1153 INCL_WARN("Returning a NKbToNKbpi channel" << '\n');
1161 isElastic = false; 1154 isElastic = false;
1162 return new NKbToNKbpiChannel( 1155 return new NKbToNKbpiChannel(particle1, particle2);
1163 } else if(quasielasticCX>0.) { 1156 } else if(quasielasticCX>0.) {
1164 INCL_WARN("Returning a quasi- 1157 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1165 isElastic = false; // true ?? 1158 isElastic = false; // true ??
1166 return new NKbToNKbChannel(pa 1159 return new NKbToNKbChannel(particle1, particle2);
1167 } else { 1160 } else {
1168 INCL_WARN("Returning an elast 1161 INCL_WARN("Returning an elastic channel" << '\n');
1169 isElastic = true; 1162 isElastic = true;
1170 return new NKbElasticChannel( 1163 return new NKbElasticChannel(particle1, particle2);
1171 } 1164 }
1172 } 1165 }
1173 } else if ((particle1->isNucleon() && par 1166 } else if ((particle1->isNucleon() && particle2->isLambda()) || (particle2->isNucleon() && particle1->isLambda())) {
1174 //// NLambda 1167 //// NLambda
1175 const G4double elasticCX = CrossSecti 1168 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1176 const G4double NLToNSCX = CrossSectio 1169 const G4double NLToNSCX = CrossSections::NLToNS(particle1,particle2);
1177 const G4double totCX=CrossSections::t 1170 const G4double totCX=CrossSections::total(particle1, particle2);
1178 // assert(std::fabs(totCX-elasticCX-NLToNSCX) 1171 // assert(std::fabs(totCX-elasticCX-NLToNSCX)<0.1);
1179 1172
1180 const G4double rChannel=Random::shoot 1173 const G4double rChannel=Random::shoot() * totCX;
1181 if(elasticCX > rChannel){ 1174 if(elasticCX > rChannel){
1182 // Elastic NLambda channel is chosen 1175 // Elastic NLambda channel is chosen
1183 isElastic = true; 1176 isElastic = true;
1184 INCL_DEBUG("NLambda interaction: 1177 INCL_DEBUG("NLambda interaction: elastic channel chosen" << '\n');
1185 return new NYElasticChannel(parti 1178 return new NYElasticChannel(particle1, particle2);
1186 } else if(elasticCX + NLToNSCX > rCha 1179 } else if(elasticCX + NLToNSCX > rChannel){
1187 // Quasi-elastic NLambda channel is chosen 1180 // Quasi-elastic NLambda channel is chosen
1188 isElastic = false; // true ?? 1181 isElastic = false; // true ??
1189 INCL_DEBUG("NLambda interaction: 1182 INCL_DEBUG("NLambda interaction: quasi-elastic channel chosen" << '\n');
1190 return new NLToNSChannel(particle 1183 return new NLToNSChannel(particle1, particle2);
1191 } else { 1184 } else {
1192 INCL_WARN("inconsistency within t 1185 INCL_WARN("inconsistency within the NLambda Cross Sections (sum!=inelastic)" << '\n');
1193 if(NLToNSCX>0.) { 1186 if(NLToNSCX>0.) {
1194 INCL_WARN("Returning a quasi- 1187 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1195 isElastic = false; // true ?? 1188 isElastic = false; // true ??
1196 return new NLToNSChannel(part 1189 return new NLToNSChannel(particle1, particle2);
1197 } else { 1190 } else {
1198 INCL_WARN("Returning an elast 1191 INCL_WARN("Returning an elastic channel" << '\n');
1199 isElastic = true; 1192 isElastic = true;
1200 return new NYElasticChannel(p 1193 return new NYElasticChannel(particle1, particle2);
1201 } 1194 }
1202 } 1195 }
1203 } else if ((particle1->isNucleon() && par 1196 } else if ((particle1->isNucleon() && particle2->isSigma()) || (particle2->isNucleon() && particle1->isSigma())) {
1204 //// NSigma 1197 //// NSigma
1205 const G4double elasticCX = CrossSecti 1198 const G4double elasticCX = CrossSections::elastic(particle1,particle2);
1206 const G4double NSToNLCX = CrossSectio 1199 const G4double NSToNLCX = CrossSections::NSToNL(particle1,particle2);
1207 const G4double NSToNSCX = CrossSectio 1200 const G4double NSToNSCX = CrossSections::NSToNS(particle1,particle2);
1208 const G4double totCX=CrossSections::t 1201 const G4double totCX=CrossSections::total(particle1, particle2);
1209 // assert(std::fabs(totCX-elasticCX-NSToNLCX- 1202 // assert(std::fabs(totCX-elasticCX-NSToNLCX-NSToNSCX)<0.1);
1210 1203
1211 const G4double rChannel=Random::shoot 1204 const G4double rChannel=Random::shoot() * totCX;
1212 if(elasticCX > rChannel){ 1205 if(elasticCX > rChannel){
1213 // Elastic NSigma channel is chosen 1206 // Elastic NSigma channel is chosen
1214 isElastic = true; 1207 isElastic = true;
1215 INCL_DEBUG("NSigma interaction: e 1208 INCL_DEBUG("NSigma interaction: elastic channel chosen" << '\n');
1216 return new NYElasticChannel(parti 1209 return new NYElasticChannel(particle1, particle2);
1217 } else if(elasticCX + NSToNLCX > rCha 1210 } else if(elasticCX + NSToNLCX > rChannel){
1218 // NSigma -> NLambda channel is chosen 1211 // NSigma -> NLambda channel is chosen
1219 isElastic = false; // true ?? 1212 isElastic = false; // true ??
1220 INCL_DEBUG("NSigma interaction: N 1213 INCL_DEBUG("NSigma interaction: NLambda channel chosen" << '\n');
1221 return new NSToNLChannel(particle 1214 return new NSToNLChannel(particle1, particle2);
1222 } else if(elasticCX + NSToNLCX + NSTo 1215 } else if(elasticCX + NSToNLCX + NSToNSCX > rChannel){
1223 // NSigma -> NSigma quasi-elastic channel is 1216 // NSigma -> NSigma quasi-elastic channel is chosen
1224 isElastic = false; // true ?? 1217 isElastic = false; // true ??
1225 INCL_DEBUG("NSigma interaction: N 1218 INCL_DEBUG("NSigma interaction: NSigma quasi-elastic channel chosen" << '\n');
1226 return new NSToNSChannel(particle 1219 return new NSToNSChannel(particle1, particle2);
1227 } else { 1220 } else {
1228 INCL_WARN("inconsistency within t 1221 INCL_WARN("inconsistency within the NSigma Cross Sections (sum!=inelastic)" << '\n');
1229 if(NSToNSCX>0.) { 1222 if(NSToNSCX>0.) {
1230 INCL_WARN("Returning a quasi- 1223 INCL_WARN("Returning a quasi-elastic channel" << '\n');
1231 isElastic = false; // true ?? 1224 isElastic = false; // true ??
1232 return new NSToNSChannel(part 1225 return new NSToNSChannel(particle1, particle2);
1233 } else if(NSToNLCX>0.) { 1226 } else if(NSToNLCX>0.) {
1234 INCL_WARN("Returning a NLambd 1227 INCL_WARN("Returning a NLambda channel" << '\n');
1235 isElastic = false; // true ?? 1228 isElastic = false; // true ??
1236 return new NSToNLChannel(part 1229 return new NSToNLChannel(particle1, particle2);
1237 } else { 1230 } else {
1238 INCL_WARN("Returning an elast 1231 INCL_WARN("Returning an elastic channel" << '\n');
1239 isElastic = true; 1232 isElastic = true;
1240 return new NYElasticChannel(p 1233 return new NYElasticChannel(particle1, particle2);
1241 } 1234 }
1242 } 1235 }
1243 } else if ((particle1->isNucleon() && par <<
1244 //// NNbar <<
1245 const G4double totCX = CrossSections: <<
1246 const G4double NNbElasticCX = CrossSe <<
1247 const G4double NNbCEXCX = CrossSectio <<
1248 const G4double NNbToLLbCX = CrossSect <<
1249 const G4double NNbToNNbpiCX = CrossSe <<
1250 const G4double NNbToNNb2piCX = CrossS <<
1251 const G4double NNbToNNb3piCX = CrossS <<
1252 const G4double AnnihilationCX = Cross <<
1253 <<
1254 // assert(std::fabs(totCX-NNbElasticCX-NNbCEX <<
1255 <<
1256 const G4double rChannel=Random::shoot <<
1257 if (NNbElasticCX > rChannel) { <<
1258 // NNbar (elastic) channel is cho <<
1259 isElastic = true; <<
1260 //INCL_WARN("NNbar interaction: N <<
1261 return new NNbarElasticChannel(pa <<
1262 } else if (NNbElasticCX + NNbCEXCX > <<
1263 // NNbar (CEX) channel is chosen <<
1264 isElastic = false; // may be char <<
1265 //INCL_WARN("NNbar interaction: N <<
1266 return new NNbarCEXChannel(partic <<
1267 } else if (NNbElasticCX + NNbCEXCX + <<
1268 // NNbarToLLbar channel is chosen <<
1269 isElastic = false; // may be char <<
1270 //INCL_WARN("NNbar interaction: N <<
1271 return new NNbarToLLbarChannel(pa <<
1272 } else if (NNbElasticCX + NNbCEXCX + <<
1273 // NNbar to NNbar pi channel is c <<
1274 isElastic = false; <<
1275 //INCL_WARN("NNbar interaction: N <<
1276 return new NNbarToNNbarpiChannel( <<
1277 } else if (NNbElasticCX + NNbCEXCX + <<
1278 // NNbar to NNbar 2pi channel is <<
1279 isElastic = false; <<
1280 //INCL_WARN("NNbar interaction: N <<
1281 return new NNbarToNNbar2piChannel <<
1282 } else if (NNbElasticCX + NNbCEXCX + <<
1283 // NNbar to NNbar 3pi channel is <<
1284 isElastic = false; <<
1285 //INCL_WARN("NNbar interaction: N <<
1286 return new NNbarToNNbar3piChannel <<
1287 } else if (NNbElasticCX + NNbCEXCX + <<
1288 // NNbar annihilation channel is <<
1289 isElastic = false; <<
1290 AnnihilationType atype; <<
1291 if((particle1->getType()==antiPro <<
1292 atype = PTypeInFlight; <<
1293 } <<
1294 else if((particle1->getType()==an <<
1295 atype = NTypeInFlight; <<
1296 } <<
1297 else if((particle1->getType()==an <<
1298 atype = NbarPTypeInFlight; <<
1299 } <<
1300 else if((particle1->getType()==an <<
1301 atype = NbarNTypeInFlight; <<
1302 } <<
1303 else{ <<
1304 atype = Def; <<
1305 INCL_ERROR("Annihilation type p <<
1306 } <<
1307 theNucleus->setAType(atype); <<
1308 return new NNbarToAnnihilationCha <<
1309 } else { <<
1310 INCL_WARN("Inconsistency within <<
1311 if (NNbToNNb3piCX > 0.0) { <<
1312 INCL_WARN("Returning an NNb <<
1313 isElastic = false; <<
1314 return new NNbarToNNbar3piC <<
1315 } else if (NNbToNNb2piCX > 0.0) <<
1316 INCL_WARN("Returning an NNb <<
1317 isElastic = false; <<
1318 return new NNbarToNNbar2piC <<
1319 } else if (NNbToNNbpiCX > 0.0) <<
1320 INCL_WARN("Returning an NNb <<
1321 isElastic = false; <<
1322 return new NNbarToNNbarpiCh <<
1323 } else if (AnnihilationCX > 0.0 <<
1324 INCL_WARN("Returning an NNb <<
1325 isElastic = false; <<
1326 AnnihilationType atype; <<
1327 if((particle1->getType()==a <<
1328 atype = PTypeInFlight; <<
1329 } <<
1330 else if((particle1->getType <<
1331 atype = NTypeInFlight; <<
1332 } <<
1333 else if((particle1->getType <<
1334 atype = NbarPTypeInFlight <<
1335 } <<
1336 else if((particle1->getType <<
1337 atype = NbarNTypeInFlight <<
1338 } <<
1339 else{ <<
1340 atype = Def; <<
1341 INCL_ERROR("Annihilation <<
1342 } <<
1343 theNucleus->setAType(atype) <<
1344 return new NNbarToAnnihilat <<
1345 } else if (NNbCEXCX > 0.0) { <<
1346 INCL_WARN("Returning an NNb <<
1347 isElastic = false; <<
1348 return new NNbarCEXChannel( <<
1349 } else if (NNbToLLbCX > 0.0) { <<
1350 INCL_WARN("Returning an NNb <<
1351 isElastic = false; <<
1352 return new NNbarToLLbarChan <<
1353 } else { <<
1354 INCL_WARN("Elastic NNbar ch <<
1355 isElastic = true; <<
1356 return new NNbarElasticChan <<
1357 } <<
1358 } <<
1359 } 1236 }
1360 << 1237
1361 else { 1238 else {
1362 INCL_DEBUG("BinaryCollisionAvatar can o 1239 INCL_DEBUG("BinaryCollisionAvatar can only handle nucleons (for the moment)."
1363 << '\n' 1240 << '\n'
1364 << particle1->print() 1241 << particle1->print()
1365 << '\n' 1242 << '\n'
1366 << particle2->print() 1243 << particle2->print()
1367 << '\n'); 1244 << '\n');
1368 InteractionAvatar::restoreParticles(); 1245 InteractionAvatar::restoreParticles();
1369 return NULL; 1246 return NULL;
1370 } 1247 }
1371 } 1248 }
1372 1249
1373 void BinaryCollisionAvatar::preInteraction( 1250 void BinaryCollisionAvatar::preInteraction() {
1374 isParticle1Spectator = particle1->isTarge 1251 isParticle1Spectator = particle1->isTargetSpectator();
1375 isParticle2Spectator = particle2->isTarge 1252 isParticle2Spectator = particle2->isTargetSpectator();
1376 InteractionAvatar::preInteraction(); 1253 InteractionAvatar::preInteraction();
1377 } 1254 }
1378 1255
1379 void BinaryCollisionAvatar::postInteraction 1256 void BinaryCollisionAvatar::postInteraction(FinalState *fs) {
1380 // Call the postInteraction method of the 1257 // Call the postInteraction method of the parent class
1381 // (provides Pauli blocking and enforces 1258 // (provides Pauli blocking and enforces energy conservation)
1382 InteractionAvatar::postInteraction(fs); 1259 InteractionAvatar::postInteraction(fs);
1383 1260
1384 switch(fs->getValidity()) { 1261 switch(fs->getValidity()) {
1385 case PauliBlockedFS: 1262 case PauliBlockedFS:
1386 theNucleus->getStore()->getBook().inc 1263 theNucleus->getStore()->getBook().incrementBlockedCollisions();
1387 break; 1264 break;
1388 case NoEnergyConservationFS: 1265 case NoEnergyConservationFS:
1389 case ParticleBelowFermiFS: 1266 case ParticleBelowFermiFS:
1390 case ParticleBelowZeroFS: 1267 case ParticleBelowZeroFS:
1391 break; 1268 break;
1392 case ValidFS: 1269 case ValidFS:
1393 Book &theBook = theNucleus->getStore( 1270 Book &theBook = theNucleus->getStore()->getBook();
1394 theBook.incrementAcceptedCollisions() 1271 theBook.incrementAcceptedCollisions();
1395 1272
1396 if(theBook.getAcceptedCollisions() == 1273 if(theBook.getAcceptedCollisions() == 1) {
1397 // Store time and cross section of 1274 // Store time and cross section of the first collision
1398 G4double t = theBook.getCurrentTime 1275 G4double t = theBook.getCurrentTime();
1399 theBook.setFirstCollisionTime(t); 1276 theBook.setFirstCollisionTime(t);
1400 theBook.setFirstCollisionXSec(oldXS 1277 theBook.setFirstCollisionXSec(oldXSec);
1401 // Increase the number of Kaon by 1 1278 // Increase the number of Kaon by 1
1402 if(isStrangeProduction) theNucleus- 1279 if(isStrangeProduction) theNucleus->setNumberOfKaon(theNucleus->getNumberOfKaon()+1);
1403 // Store position and momentum of t 1280 // Store position and momentum of the spectator on the first
1404 // collision 1281 // collision
1405 if((isParticle1Spectator && isParti 1282 if((isParticle1Spectator && isParticle2Spectator) || (!isParticle1Spectator && !isParticle2Spectator)) {
1406 INCL_ERROR("First collision must 1283 INCL_ERROR("First collision must be within a target spectator and a non-target spectator");
1407 } 1284 }
1408 if(isParticle1Spectator) { 1285 if(isParticle1Spectator) {
1409 theBook.setFirstCollisionSpectato 1286 theBook.setFirstCollisionSpectatorPosition(backupParticle1->getPosition().mag());
1410 theBook.setFirstCollisionSpectato 1287 theBook.setFirstCollisionSpectatorMomentum(backupParticle1->getMomentum().mag());
1411 } else { 1288 } else {
1412 theBook.setFirstCollisionSpectato 1289 theBook.setFirstCollisionSpectatorPosition(backupParticle2->getPosition().mag());
1413 theBook.setFirstCollisionSpectato 1290 theBook.setFirstCollisionSpectatorMomentum(backupParticle2->getMomentum().mag());
1414 } 1291 }
1415 1292
1416 // Store the elasticity of the firs 1293 // Store the elasticity of the first collision
1417 theBook.setFirstCollisionIsElastic( 1294 theBook.setFirstCollisionIsElastic(isElastic);
1418 } 1295 }
1419 } 1296 }
1420 return; 1297 return;
1421 } 1298 }
1422 1299
1423 std::string BinaryCollisionAvatar::dump() c 1300 std::string BinaryCollisionAvatar::dump() const {
1424 std::stringstream ss; 1301 std::stringstream ss;
1425 ss << "(avatar " << theTime <<" 'nn-colli 1302 ss << "(avatar " << theTime <<" 'nn-collision" << '\n'
1426 << "(list " << '\n' 1303 << "(list " << '\n'
1427 << particle1->dump() 1304 << particle1->dump()
1428 << particle2->dump() 1305 << particle2->dump()
1429 << "))" << '\n'; 1306 << "))" << '\n';
1430 return ss.str(); 1307 return ss.str();
1431 } 1308 }
1432 1309
1433 } 1310 }
1434 1311