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