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