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 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics 28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA 29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA 31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege 32 // Davide Mancusi, CEA-Saclay, France << 33 // 32 // 34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1 35 34 36 #include "globals.hh" 35 #include "globals.hh" 37 36 38 /* 37 /* 39 * G4INCLBinaryCollisionAvatar.cc 38 * G4INCLBinaryCollisionAvatar.cc 40 * 39 * 41 * \date Jun 5, 2009 40 * \date Jun 5, 2009 42 * \author Pekka Kaitaniemi 41 * \author Pekka Kaitaniemi 43 */ 42 */ 44 43 45 #include "G4INCLBinaryCollisionAvatar.hh" 44 #include "G4INCLBinaryCollisionAvatar.hh" 46 #include "G4INCLElasticChannel.hh" 45 #include "G4INCLElasticChannel.hh" 47 #include "G4INCLRecombinationChannel.hh" 46 #include "G4INCLRecombinationChannel.hh" 48 #include "G4INCLDeltaProductionChannel.hh" 47 #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.h << 55 #include "G4INCLNNOmegaToMultiPionsChannel.hh" << 56 #include "G4INCLCrossSections.hh" 48 #include "G4INCLCrossSections.hh" 57 #include "G4INCLKinematicsUtils.hh" 49 #include "G4INCLKinematicsUtils.hh" 58 #include "G4INCLRandom.hh" 50 #include "G4INCLRandom.hh" 59 #include "G4INCLParticleTable.hh" 51 #include "G4INCLParticleTable.hh" 60 #include "G4INCLPauliBlocking.hh" 52 #include "G4INCLPauliBlocking.hh" 61 #include "G4INCLPiNElasticChannel.hh" << 53 #include "G4INCLPionNucleonChannel.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. << 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" << 109 #include "G4INCLStore.hh" 54 #include "G4INCLStore.hh" 110 #include "G4INCLBook.hh" 55 #include "G4INCLBook.hh" 111 #include "G4INCLLogger.hh" 56 #include "G4INCLLogger.hh" 112 #include <string> 57 #include <string> 113 #include <sstream> 58 #include <sstream> 114 // #include <cassert> 59 // #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 60 123 namespace G4INCL { 61 namespace G4INCL { 124 62 125 // WARNING: if you update the default cutNN << 63 const G4double BinaryCollisionAvatar::cutNN = 1910; 126 // cutNNSquared variable, too. << 64 const G4double BinaryCollisionAvatar::cutNNSquared = cutNN*cutNN; 127 G4ThreadLocal G4double BinaryCollisionAvatar << 128 G4ThreadLocal G4double BinaryCollisionAvatar << 129 G4ThreadLocal G4double BinaryCollisionAvatar << 130 65 131 BinaryCollisionAvatar::BinaryCollisionAvatar 66 BinaryCollisionAvatar::BinaryCollisionAvatar(G4double time, G4double crossSection, 132 G4INCL::Nucleus *n, G4INCL::Particle *p1 67 G4INCL::Nucleus *n, G4INCL::Particle *p1, G4INCL::Particle *p2) 133 : InteractionAvatar(time, n, p1, p2), theC << 68 : InteractionAvatar(time, n, p1, p2), theCrossSection(crossSection) 134 isParticle1Spectator(false), << 135 isParticle2Spectator(false), << 136 isElastic(false), << 137 isStrangeProduction(false) << 138 { 69 { 139 setType(CollisionAvatarType); 70 setType(CollisionAvatarType); 140 } 71 } 141 72 142 BinaryCollisionAvatar::~BinaryCollisionAvata 73 BinaryCollisionAvatar::~BinaryCollisionAvatar() { 143 } 74 } 144 75 145 G4INCL::IChannel* BinaryCollisionAvatar::get << 76 G4INCL::IChannel* BinaryCollisionAvatar::getChannel() const { 146 // We already check cutNN at avatar creati 77 // We already check cutNN at avatar creation time, but we have to check it 147 // again here. For composite projectiles, 78 // again here. For composite projectiles, we might have created independent 148 // avatars with no cutNN before any collis 79 // avatars with no cutNN before any collision took place. 149 if(particle1->isNucleon() 80 if(particle1->isNucleon() 150 && particle2->isNucleon() 81 && particle2->isNucleon() 151 && theNucleus->getStore()->getBook().g << 82 && theNucleus->getStore()->getBook()->getAcceptedCollisions()!=0) { 152 const G4double energyCM2 = KinematicsUti 83 const G4double energyCM2 = KinematicsUtils::squareTotalEnergyInCM(particle1, particle2); 153 // Below a certain cut value we don't do 84 // Below a certain cut value we don't do anything: 154 if(energyCM2 < cutNNSquared) { 85 if(energyCM2 < cutNNSquared) { 155 INCL_DEBUG("CM energy = sqrt(" << ener << 86 DEBUG("CM energy = sqrt(" << energyCM2 << ") MeV < sqrt(" << cutNNSquared 156 << ") MeV = cutNN" << "; returning << 87 << ") MeV = cutNN" << "; returning a NULL channel" << std::endl); 157 InteractionAvatar::restoreParticles(); 88 InteractionAvatar::restoreParticles(); 158 return NULL; 89 return NULL; 159 } 90 } 160 } 91 } 161 92 162 /** Check again the distance of approach. 93 /** Check again the distance of approach. In order for the avatar to be 163 * realised, we have to perform a check in 94 * realised, we have to perform a check in the CM system. We define a 164 * distance four-vector as 95 * distance four-vector as 165 * \f[ (0, \Delta\vec{x}), \f] 96 * \f[ (0, \Delta\vec{x}), \f] 166 * where \f$\Delta\vec{x}\f$ is the distan 97 * where \f$\Delta\vec{x}\f$ is the distance vector of the particles at 167 * their minimum distance of approach (i.e 98 * their minimum distance of approach (i.e. at the avatar time). By 168 * boosting this four-vector to the CM fra 99 * boosting this four-vector to the CM frame of the two particles and we 169 * obtain a new four vector 100 * obtain a new four vector 170 * \f[ (\Delta t', \Delta\vec{x}'), \f] 101 * \f[ (\Delta t', \Delta\vec{x}'), \f] 171 * with a non-zero time component (the col 102 * with a non-zero time component (the collision happens simultaneously for 172 * the two particles in the lab system, bu 103 * the two particles in the lab system, but not in the CM system). In order 173 * for the avatar to be realised, we requi 104 * for the avatar to be realised, we require that 174 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/ 105 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/\pi}.\f] 175 * Note that \f$|\Delta\vec{x}'|\leq|\Delt 106 * Note that \f$|\Delta\vec{x}'|\leq|\Delta\vec{x}|\f$; thus, the condition 176 * above is more restrictive than the chec 107 * above is more restrictive than the check that we perform in 177 * G4INCL::Propagation::StandardPropagatio 108 * G4INCL::Propagation::StandardPropagationModel::generateBinaryCollisionAvatar. 178 * In other words, the avatar generation c 109 * In other words, the avatar generation cannot miss any physical collision 179 * avatars. 110 * avatars. 180 */ 111 */ 181 ThreeVector minimumDistance = particle1->g 112 ThreeVector minimumDistance = particle1->getPosition(); 182 minimumDistance -= particle2->getPosition( 113 minimumDistance -= particle2->getPosition(); 183 const G4double betaDotX = boostVector.dot( 114 const G4double betaDotX = boostVector.dot(minimumDistance); 184 const G4double minDist = Math::tenPi*(mini 115 const G4double minDist = Math::tenPi*(minimumDistance.mag2() + betaDotX*betaDotX / (1.-boostVector.mag2())); 185 if(minDist > theCrossSection) { 116 if(minDist > theCrossSection) { 186 INCL_DEBUG("CM distance of approach is t << 117 DEBUG("CM distance of approach is too small: " << minDist << ">" << 187 theCrossSection <<"; returning a NULL << 118 theCrossSection <<"; returning a NULL channel" << std::endl); 188 InteractionAvatar::restoreParticles(); 119 InteractionAvatar::restoreParticles(); 189 return NULL; 120 return NULL; 190 } 121 } 191 122 192 /** Bias apply for this reaction in order << 123 if(particle1->isNucleon() && particle2->isNucleon()) { // NN->NN 193 * ParticleBias for all stange particles. << 124 G4double elasticCX = CrossSections::elastic(particle1, 194 * Can be reduced after because of the saf << 125 particle2); 195 */ << 126 G4double deltaProductionCX = CrossSections::deltaProduction(particle1, 196 G4double bias_apply = 1.; << 127 particle2); 197 if(bias != 1.) bias_apply = Particle::getB << 128 198 << 129 G4bool isElastic = true; 199 //// NN << 130 if(elasticCX/(elasticCX + deltaProductionCX) < Random::shoot()) { 200 if(particle1->isNucleon() && particle2->is << 131 // NN -> N Delta channel is chosen 201 << 202 G4double NLKProductionCX = CrossSections << 203 G4double NSKProductionCX = CrossSections << 204 G4double NLKpiProductionCX = CrossSectio << 205 G4double NSKpiProductionCX = CrossSectio << 206 G4double NLK2piProductionCX = CrossSecti << 207 G4double NSK2piProductionCX = CrossSecti << 208 G4double NNKKbProductionCX = CrossSectio << 209 G4double NNMissingCX = CrossSections::NN << 210 << 211 const G4double UnStrangeProdCX = CrossSe << 212 + CrossSect << 213 + CrossSect << 214 + CrossSect << 215 + CrossSect << 216 + CrossSect << 217 const G4double StrangenessProdCX = (NLKP << 218 << 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 << 256 << 257 // assert(std::fabs(totCX-elasticCX-deltaProdu << 258 << 259 const G4double rChannel=Random::shoot() << 260 << 261 if(elasticCX > rChannel) { << 262 // Elastic NN channel << 263 isElastic = true; << 264 INCL_DEBUG("NN interaction: elastic ch << 265 weight = counterweight; << 266 return new ElasticChannel(particle1, p << 267 } else if((elasticCX + deltaProductionCX << 268 isElastic = false; << 269 // NN -> N Delta channel is chosen << 270 INCL_DEBUG("NN interaction: Delta chan << 271 weight = counterweight; << 272 return new DeltaProductionChannel(part << 273 } else if(elasticCX + deltaProductionCX << 274 isElastic = false; << 275 // NN -> PiNN channel is chosen << 276 INCL_DEBUG("NN interaction: one Pion c << 277 weight = counterweight; << 278 return new NNToMultiPionsChannel(1,par << 279 } else if(elasticCX + deltaProductionCX << 280 isElastic = false; << 281 // NN -> 2PiNN channel is chosen << 282 INCL_DEBUG("NN interaction: two Pions << 283 weight = counterweight; << 284 return new NNToMultiPionsChannel(2,par << 285 } else if(elasticCX + deltaProductionCX << 286 isElastic = false; << 287 // NN -> 3PiNN channel is chosen << 288 INCL_DEBUG("NN interaction: three Pion << 289 weight = counterweight; << 290 return new NNToMultiPionsChannel(3,par << 291 } else if(elasticCX + deltaProductionCX << 292 isElastic = false; << 293 // NN -> 4PiNN channel is chosen << 294 INCL_DEBUG("NN interaction: four << 295 weight = counterweight; << 296 return new NNToMultiPionsChannel << 297 } else if(elasticCX + deltaProductionCX << 298 + etaProductionCX > << 299 isElastic = false; << 300 // NN -> NNEta channel is chosen << 301 INCL_DEBUG("NN interaction: Eta channel << 302 weight = counterweight; << 303 return new NNToNNEtaChannel(particle1, << 304 } else if(elasticCX + deltaProductionCX << 305 + etaProductionCX + << 306 isElastic = false; << 307 // NN -> N Delta Eta channel is chosen << 308 INCL_DEBUG("NN interaction: Delta Eta c << 309 weight = counterweight; << 310 return new NDeltaEtaProductionChannel(p << 311 } else if(elasticCX + deltaProductionCX << 312 + etaProductionCX + << 313 isElastic = false; << 314 // NN -> EtaPiNN channel is chosen << 315 INCL_DEBUG("NN interaction: Eta + one P << 316 weight = counterweight; << 317 return new NNEtaToMultiPionsChannel(1,p << 318 } else if(elasticCX + deltaProductionCX << 319 + etaProductionCX + << 320 isElastic = false; << 321 // NN -> Eta2PiNN channel is chosen << 322 INCL_DEBUG("NN interaction: Eta + two P << 323 weight = counterweight; << 324 return new NNEtaToMultiPionsChannel(2,p << 325 } else if(elasticCX + deltaProductionCX << 326 + etaProductionCX + << 327 isElastic = false; << 328 // NN -> Eta3PiNN channel is chosen << 329 INCL_DEBUG("NN interaction: Eta + three << 330 weight = counterweight; << 331 return new NNEtaToMultiPionsChannel(3,p << 332 } else if(elasticCX + deltaProductionCX << 333 + etaProductionCX + << 334 isElastic = false; << 335 // NN -> Eta4PiNN channel is chosen << 336 INCL_DEBUG("NN interaction: Eta + four << 337 weight = counterweight; << 338 return new NNEtaToMultiPionsChannel(4,p << 339 } else if(elasticCX + deltaProductionCX << 340 + etaProductionCX + << 341 + omegaProductionCX << 342 isElastic = false; << 343 // NN -> NNOmega channel is chosen << 344 INCL_DEBUG("NN interaction: Omega chann << 345 weight = counterweight; << 346 return new NNToNNOmegaChannel(particle1 << 347 } else if(elasticCX + deltaProductionCX << 348 + etaProductionCX + << 349 + omegaProductionCX << 350 isElastic = false; << 351 // NN -> N Delta Omega channel is chosen << 352 INCL_DEBUG("NN interaction: Delta Omega << 353 weight = counterweight; << 354 return new NDeltaOmegaProductionChannel << 355 } else if(elasticCX + deltaProductionCX << 356 + etaProductionCX + << 357 + omegaProductionCX << 358 isElastic = false; << 359 // NN -> OmegaPiNN channel is chosen << 360 INCL_DEBUG("NN interaction: Omega + one << 361 weight = counterweight; << 362 return new NNOmegaToMultiPionsChannel(1 << 363 } else if(elasticCX + deltaProductionCX << 364 + etaProductionCX + << 365 + omegaProductionCX << 366 isElastic = false; << 367 // NN -> Omega2PiNN channel is chosen << 368 INCL_DEBUG("NN interaction: Omega + two << 369 weight = counterweight; << 370 return new NNOmegaToMultiPionsChannel(2 << 371 } else if(elasticCX + deltaProductionCX << 372 + etaProductionCX + << 373 + omegaProductionCX << 374 isElastic = false; << 375 // NN -> Omega3PiNN channel is chosen << 376 INCL_DEBUG("NN interaction: Omega + thr << 377 weight = counterweight; << 378 return new NNOmegaToMultiPionsChannel(3 << 379 } else if(elasticCX + deltaProductionCX << 380 + etaProductionCX + << 381 + omegaProductionCX << 382 isElastic = false; << 383 // NN -> Omega4PiNN channel is chosen << 384 INCL_DEBUG("NN interaction: Omega + fou << 385 weight = counterweight; << 386 return new NNOmegaToMultiPionsChannel(4 << 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; 132 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 { << 468 INCL_WARN("inconsistency within the NN << 469 if(NNMissingCX>0.) { << 470 INCL_WARN("Returning an Missing St << 471 weight = bias_apply; << 472 isElastic = false; << 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.) << 523 INCL_WARN("Returning an Omega + th << 524 weight = counterweight; << 525 isElastic = false; << 526 return new NNOmegaToMultiPionsChan << 527 } else if(omegatwoPiProductionCX>0.) { << 528 INCL_WARN("Returning an Omega + tw << 529 weight = counterweight; << 530 isElastic = false; << 531 return new NNOmegaToMultiPionsChan << 532 } else if(omegaonePiProductionCX>0.) { << 533 INCL_WARN("Returning an Omega + on << 534 weight = counterweight; << 535 isElastic = false; << 536 return new NNOmegaToMultiPionsChannel << 537 } else if(omegadeltaProductionCX>0.) { << 538 INCL_WARN("Returning an Omega + De << 539 weight = counterweight; << 540 isElastic = false; << 541 return new NDeltaOmegaProductionCh << 542 } else if(omegaProductionCX>0.) { << 543 INCL_WARN("Returning an Omega chan << 544 weight = counterweight; << 545 isElastic = false; << 546 return new NNToNNOmegaChannel(part << 547 } else if(etafourPiProductionCX>0.) { << 548 INCL_WARN("Returning an Eta + four << 549 weight = counterweight; << 550 isElastic = false; << 551 return new NNEtaToMultiPionsChanne << 552 } else if(etathreePiProductionCX>0.) { << 553 INCL_WARN("Returning an Eta + thre << 554 weight = counterweight; << 555 isElastic = false; << 556 return new NNEtaToMultiPionsChanne << 557 } else if(etatwoPiProductionCX>0.) { << 558 INCL_WARN("Returning an Eta + two << 559 weight = counterweight; << 560 isElastic = false; << 561 return new NNEtaToMultiPionsChanne << 562 } else if(etaonePiProductionCX>0.) { << 563 INCL_WARN("Returning an Eta + one << 564 weight = counterweight; << 565 isElastic = false; << 566 return new NNEtaToMultiPionsChanne << 567 } else if(etadeltaProductionCX>0.) { << 568 INCL_WARN("Returning an Eta + Delt << 569 weight = counterweight; << 570 isElastic = false; << 571 return new NDeltaEtaProductionChan << 572 } else if(etaProductionCX>0.) { << 573 INCL_WARN("Returning an Eta channe << 574 weight = counterweight; << 575 isElastic = false; << 576 return new NNToNNEtaChannel(partic << 577 } else if(fourPiProductionCX>0.) { << 578 INCL_WARN("Returning a 4pi channel << 579 weight = counterweight; << 580 isElastic = false; << 581 return new NNToMultiPionsChannel(4 << 582 } else if(threePiProductionCX>0.) { << 583 INCL_WARN("Returning a 3pi channel << 584 weight = counterweight; << 585 isElastic = false; << 586 return new NNToMultiPionsChannel(3 << 587 } else if(twoPiProductionCX>0.) { << 588 INCL_WARN("Returning a 2pi channel << 589 weight = counterweight; << 590 isElastic = false; << 591 return new NNToMultiPionsChannel(2 << 592 } else if(onePiProductionCX>0.) { << 593 INCL_WARN("Returning a 1pi channel << 594 weight = counterweight; << 595 isElastic = false; << 596 return new NNToMultiPionsChannel(1 << 597 } else if(deltaProductionCX>0.) { << 598 INCL_WARN("Returning a delta-produ << 599 weight = counterweight; << 600 isElastic = false; << 601 return new DeltaProductionChannel( << 602 } else { << 603 INCL_WARN("Returning an elastic ch << 604 weight = counterweight; << 605 isElastic = true; << 606 return new ElasticChannel(particle << 607 } << 608 } << 609 << 610 //// NDelta << 611 } << 612 else if((particle1->isNucleon() && particl << 613 (particle1->isDelta() && part << 614 << 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 << 642 if(elasticCX > rChannel) { << 643 isElastic = true; << 644 // Elastic N Delta channel << 645 INCL_DEBUG("NDelta interaction: e << 646 weight = counterweight; << 647 return new ElasticChannel(particl << 648 } else if (elasticCX + recombination << 649 isElastic = false; << 650 // Recombination << 651 // NDelta -> NN channel is chosen << 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 << 698 //// DeltaDelta << 699 } else if(particle1->isDelta() && particle << 700 isElastic = true; << 701 INCL_DEBUG("DeltaDelta interaction: el << 702 return new ElasticChannel(particle1, p << 703 << 704 //// PiN << 705 } else if(isPiN) { << 706 << 707 G4double LKProdCX = CrossSections::NpiTo << 708 G4double SKProdCX = CrossSections::NpiTo << 709 G4double LKpiProdCX = CrossSections::Npi << 710 G4double SKpiProdCX = CrossSections::Npi << 711 G4double LK2piProdCX = CrossSections::Np << 712 G4double SK2piProdCX = CrossSections::Np << 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 } 133 } 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 << 746 << 747 // assert(std::fabs(totCX-elasticCX-deltaProdu << 748 134 749 const G4double rChannel=Random::shoot() << 135 if(isElastic) { // Elastic NN channel 750 << 136 DEBUG("NN interaction: elastic channel chosen" << std::endl); 751 if(elasticCX > rChannel) { << 137 return new ElasticChannel(theNucleus, particle1, particle2); 752 isElastic = true; << 138 } else { // Delta production 753 // Elastic PiN channel << 139 // Inelastic NN channel 754 INCL_DEBUG("PiN interaction: elastic c << 140 DEBUG("NN interaction: inelastic channel chosen" << std::endl); 755 weight = counterweight; << 141 return new DeltaProductionChannel(particle1, particle2, theNucleus); 756 return new PiNElasticChannel(particle1 << 142 } 757 } else if(elasticCX + deltaProductionCX << 143 } else if((particle1->isNucleon() && particle2->isDelta()) || 758 isElastic = false; << 144 (particle1->isDelta() && particle2->isNucleon())) { 759 // PiN -> Delta channel is chosen << 145 G4double elasticCX = CrossSections::elastic(particle1, 760 INCL_DEBUG("PiN interaction: Delta cha << 146 particle2); 761 weight = counterweight; << 147 G4double recombinationCX = CrossSections::recombination(particle1, 762 return new PiNToDeltaChannel(particle1 << 148 particle2); 763 } else if(elasticCX + deltaProductionCX << 149 764 isElastic = false; << 150 G4bool isElastic = true; 765 // PiN -> PiNPi channel is chosen << 151 if(elasticCX/(elasticCX + recombinationCX) < Random::shoot()) { 766 INCL_DEBUG("PiN interaction: one Pion << 152 // N Delta -> NN channel is chosen 767 weight = counterweight; << 768 return new PiNToMultiPionsChannel(2,pa << 769 } else if(elasticCX + deltaProductionCX << 770 isElastic = false; << 771 // PiN -> PiN2Pi channel is chosen << 772 INCL_DEBUG("PiN interaction: two Pions << 773 weight = counterweight; << 774 return new PiNToMultiPionsChannel(3,pa << 775 } else if(elasticCX + deltaProductionCX << 776 isElastic = false; << 777 // PiN -> PiN3Pi channel is chosen << 778 INCL_DEBUG("PiN interaction: three Pio << 779 weight = counterweight; << 780 return new PiNToMultiPionsChannel(4,pa << 781 } else if(elasticCX + deltaProductionCX << 782 isElastic = false; << 783 // PiN -> EtaN channel is chosen << 784 INCL_DEBUG("PiN interaction: Eta chann << 785 weight = counterweight; << 786 return new PiNToEtaChannel(particle1, << 787 } else if(elasticCX + deltaProductionCX << 788 isElastic = false; << 789 // PiN -> OmegaN channel is chosen << 790 INCL_DEBUG("PiN interaction: Omega cha << 791 weight = counterweight; << 792 return new PiNToOmegaChannel(particle1 << 793 } else if(elasticCX + deltaProductionCX << 794 + LKProdCX > rChanne << 795 isElastic = false; << 796 isStrangeProduction = true; << 797 // PiN -> LK channel is chosen << 798 INCL_DEBUG("PiN interaction: LK channe << 799 weight = bias_apply; << 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; 153 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.) { << 914 INCL_WARN("Returning a Eta channel << 915 weight = counterweight; << 916 isElastic = false; << 917 return new PiNToEtaChannel(particl << 918 } else if(threePiProductionCX>0.) { << 919 INCL_WARN("Returning a 3pi channel << 920 weight = counterweight; << 921 isElastic = false; << 922 return new PiNToMultiPionsChannel( << 923 } else if(twoPiProductionCX>0.) { << 924 INCL_WARN("Returning a 2pi channel << 925 weight = counterweight; << 926 isElastic = false; << 927 return new PiNToMultiPionsChannel( << 928 } else if(onePiProductionCX>0.) { << 929 INCL_WARN("Returning a 1pi channel << 930 weight = counterweight; << 931 isElastic = false; << 932 return new PiNToMultiPionsChannel( << 933 } else if(deltaProductionCX>0.) { << 934 INCL_WARN("Returning a delta-produ << 935 weight = counterweight; << 936 isElastic = false; << 937 return new PiNToDeltaChannel(parti << 938 } else { << 939 INCL_WARN("Returning an elastic ch << 940 weight = counterweight; << 941 isElastic = true; << 942 return new PiNElasticChannel(parti << 943 } << 944 } 154 } 945 } else if ((particle1->isNucleon() && part << 946 //// EtaN << 947 << 948 const G4double elasticCX = << 949 const G4double onePiProduc << 950 const G4double twoPiProduc << 951 const G4double totCX=Cross << 952 // assert(std::fabs(totCX-elasticCX-onePiProdu << 953 << 954 const G4double rChannel=Ra << 955 << 956 if(elasticCX > rChannel) { << 957 // Elastic EtaN channel << 958 isElastic = true; << 959 INCL_DEBUG("EtaN inter << 960 return new EtaNElastic << 961 } else if(elasticCX + oneP << 962 isElastic = false; << 963 // EtaN -> EtaPiN channel is chosen << 964 INCL_DEBUG("EtaN inter << 965 return new EtaNToPiNCh << 966 } else if(elasticCX + oneP << 967 isElastic = false; << 968 // EtaN -> EtaPiPiN channel is chosen << 969 INCL_DEBUG("EtaN inter << 970 return new EtaNToPiPiN << 971 } << 972 155 973 else { << 156 if(isElastic) { // Elastic N Delta channel 974 INCL_WARN("inconsisten << 157 DEBUG("NDelta interaction: elastic channel chosen" << std::endl); 975 if(twoPiProductionCX>0 << 158 return new ElasticChannel(theNucleus, particle1, particle2); 976 INCL_WARN("Returni << 159 } else { // Recombination 977 isElastic = false; << 160 DEBUG("NDelta interaction: recombination channel chosen" << std::endl); 978 return new EtaNToP << 161 return new RecombinationChannel(theNucleus, particle1, particle2); 979 } else if(onePiProduct << 162 } 980 INCL_WARN("Returni << 163 } else if(particle1->isDelta() && particle2->isDelta()) { 981 isElastic = false; << 164 DEBUG("DeltaDelta interaction: elastic channel chosen" << std::endl); 982 return new EtaNToP << 165 return new ElasticChannel(theNucleus, particle1, particle2); 983 } else { << 166 } else if((particle1->isNucleon() && particle2->isPion()) || 984 INCL_WARN("Returni << 167 (particle1->isPion() && particle2->isNucleon())) { 985 isElastic = true; << 168 return new PionNucleonChannel(particle1, particle2, theNucleus, shouldUseLocalEnergy()); 986 return new EtaNEla << 169 } else { 987 } << 170 DEBUG("BinaryCollisionAvatar can only handle nucleons (for the moment)." 988 } << 171 << std::endl 989 << 172 << particle1->print() 990 } else if ((particle1->isNucleon() && part << 173 << std::endl 991 //// OmegaN << 174 << particle2->print() 992 << 175 << std::endl); 993 const G4double elasticCX = << 994 const G4double onePiProduc << 995 const G4double twoPiProduc << 996 const G4double totCX=Cross << 997 // assert(std::fabs(totCX-elasticCX-onePiProdu << 998 << 999 const G4double rChannel=Ra << 1000 << 1001 if(elasticCX > rChannel) << 1002 // Elastic OmegaN channel << 1003 isElastic = true; << 1004 INCL_DEBUG("OmegaN in << 1005 return new OmegaNElas << 1006 } else if(elasticCX + one << 1007 isElastic = false; << 1008 // OmegaN -> PiN channel is chosen << 1009 INCL_DEBUG("OmegaN interaction: P << 1010 return new OmegaNToPiNChannel(par << 1011 } else if(elasticCX + one << 1012 isElastic = false; << 1013 // OmegaN -> PiPiN channel is chosen << 1014 INCL_DEBUG("OmegaN in << 1015 return new OmegaNToPi << 1016 } << 1017 else { << 1018 INCL_WARN("inconsiste << 1019 if(twoPiProductionCX> << 1020 INCL_WARN("Return << 1021 isElastic = false << 1022 return new OmegaN << 1023 } else if(onePiProduc << 1024 INCL_WARN("Return << 1025 isElastic = false << 1026 return new OmegaN << 1027 } else { << 1028 INCL_WARN("Return << 1029 isElastic = true; << 1030 return new OmegaN << 1031 } << 1032 } << 1033 } else if ((particle1->isNucleon() && par << 1034 //// KN << 1035 const G4double elasticCX = CrossSecti << 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 << 1363 << '\n' << 1364 << particle1->print() << 1365 << '\n' << 1366 << particle2->print() << 1367 << '\n'); << 1368 InteractionAvatar::restoreParticles(); 176 InteractionAvatar::restoreParticles(); 1369 return NULL; 177 return NULL; 1370 } 178 } 1371 } 179 } 1372 180 1373 void BinaryCollisionAvatar::preInteraction( 181 void BinaryCollisionAvatar::preInteraction() { 1374 isParticle1Spectator = particle1->isTarge << 1375 isParticle2Spectator = particle2->isTarge << 1376 InteractionAvatar::preInteraction(); 182 InteractionAvatar::preInteraction(); 1377 } 183 } 1378 184 1379 void BinaryCollisionAvatar::postInteraction << 185 FinalState *BinaryCollisionAvatar::postInteraction(FinalState *fs) { 1380 // Call the postInteraction method of the 186 // Call the postInteraction method of the parent class 1381 // (provides Pauli blocking and enforces 187 // (provides Pauli blocking and enforces energy conservation) 1382 InteractionAvatar::postInteraction(fs); << 188 fs = InteractionAvatar::postInteraction(fs); 1383 189 1384 switch(fs->getValidity()) { 190 switch(fs->getValidity()) { 1385 case PauliBlockedFS: 191 case PauliBlockedFS: 1386 theNucleus->getStore()->getBook().inc << 192 theNucleus->getStore()->getBook()->incrementBlockedCollisions(); 1387 break; 193 break; 1388 case NoEnergyConservationFS: 194 case NoEnergyConservationFS: 1389 case ParticleBelowFermiFS: 195 case ParticleBelowFermiFS: 1390 case ParticleBelowZeroFS: 196 case ParticleBelowZeroFS: 1391 break; 197 break; 1392 case ValidFS: 198 case ValidFS: 1393 Book &theBook = theNucleus->getStore( << 199 theNucleus->getStore()->getBook()->incrementAcceptedCollisions(); 1394 theBook.incrementAcceptedCollisions() << 200 if(theNucleus->getStore()->getBook()->getAcceptedCollisions() == 1) { 1395 << 201 G4double t = theNucleus->getStore()->getBook()->getCurrentTime(); 1396 if(theBook.getAcceptedCollisions() == << 202 theNucleus->getStore()->getBook()->setFirstCollisionTime(t); 1397 // Store time and cross section of << 203 theNucleus->getStore()->getBook()->setFirstCollisionXSec(oldXSec); 1398 G4double t = theBook.getCurrentTime << 1399 theBook.setFirstCollisionTime(t); << 1400 theBook.setFirstCollisionXSec(oldXS << 1401 // Increase the number of Kaon by 1 << 1402 if(isStrangeProduction) theNucleus- << 1403 // Store position and momentum of t << 1404 // collision << 1405 if((isParticle1Spectator && isParti << 1406 INCL_ERROR("First collision must << 1407 } << 1408 if(isParticle1Spectator) { << 1409 theBook.setFirstCollisionSpectato << 1410 theBook.setFirstCollisionSpectato << 1411 } else { << 1412 theBook.setFirstCollisionSpectato << 1413 theBook.setFirstCollisionSpectato << 1414 } << 1415 << 1416 // Store the elasticity of the firs << 1417 theBook.setFirstCollisionIsElastic( << 1418 } 204 } 1419 } 205 } 1420 return; << 206 return fs; 1421 } 207 } 1422 208 1423 std::string BinaryCollisionAvatar::dump() c 209 std::string BinaryCollisionAvatar::dump() const { 1424 std::stringstream ss; 210 std::stringstream ss; 1425 ss << "(avatar " << theTime <<" 'nn-colli << 211 ss << "(avatar " << theTime <<" 'nn-collision" << std::endl 1426 << "(list " << '\n' << 212 << "(list " << std::endl 1427 << particle1->dump() 213 << particle1->dump() 1428 << particle2->dump() 214 << particle2->dump() 1429 << "))" << '\n'; << 215 << "))" << std::endl; 1430 return ss.str(); 216 return ss.str(); 1431 } 217 } 1432 218 1433 } 219 } 1434 220