Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 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 H 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.h 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. 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" 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 126 // cutNNSquared variable, too. 127 G4ThreadLocal G4double BinaryCollisionAvatar 128 G4ThreadLocal G4double BinaryCollisionAvatar 129 G4ThreadLocal G4double BinaryCollisionAvatar 130 131 BinaryCollisionAvatar::BinaryCollisionAvatar 132 G4INCL::Nucleus *n, G4INCL::Particle *p1 133 : InteractionAvatar(time, n, p1, p2), theC 134 isParticle1Spectator(false), 135 isParticle2Spectator(false), 136 isElastic(false), 137 isStrangeProduction(false) 138 { 139 setType(CollisionAvatarType); 140 } 141 142 BinaryCollisionAvatar::~BinaryCollisionAvata 143 } 144 145 G4INCL::IChannel* BinaryCollisionAvatar::get 146 // We already check cutNN at avatar creati 147 // again here. For composite projectiles, 148 // avatars with no cutNN before any collis 149 if(particle1->isNucleon() 150 && particle2->isNucleon() 151 && theNucleus->getStore()->getBook().g 152 const G4double energyCM2 = KinematicsUti 153 // Below a certain cut value we don't do 154 if(energyCM2 < cutNNSquared) { 155 INCL_DEBUG("CM energy = sqrt(" << ener 156 << ") MeV = cutNN" << "; returning 157 InteractionAvatar::restoreParticles(); 158 return NULL; 159 } 160 } 161 162 /** Check again the distance of approach. 163 * realised, we have to perform a check in 164 * distance four-vector as 165 * \f[ (0, \Delta\vec{x}), \f] 166 * where \f$\Delta\vec{x}\f$ is the distan 167 * their minimum distance of approach (i.e 168 * boosting this four-vector to the CM fra 169 * obtain a new four vector 170 * \f[ (\Delta t', \Delta\vec{x}'), \f] 171 * with a non-zero time component (the col 172 * the two particles in the lab system, bu 173 * for the avatar to be realised, we requi 174 * \f[ |\Delta\vec{x}'| \leq \sqrt{\sigma/ 175 * Note that \f$|\Delta\vec{x}'|\leq|\Delt 176 * above is more restrictive than the chec 177 * G4INCL::Propagation::StandardPropagatio 178 * In other words, the avatar generation c 179 * avatars. 180 */ 181 ThreeVector minimumDistance = particle1->g 182 minimumDistance -= particle2->getPosition( 183 const G4double betaDotX = boostVector.dot( 184 const G4double minDist = Math::tenPi*(mini 185 if(minDist > theCrossSection) { 186 INCL_DEBUG("CM distance of approach is t 187 theCrossSection <<"; returning a NULL 188 InteractionAvatar::restoreParticles(); 189 return NULL; 190 } 191 192 /** Bias apply for this reaction in order 193 * ParticleBias for all stange particles. 194 * Can be reduced after because of the saf 195 */ 196 G4double bias_apply = 1.; 197 if(bias != 1.) bias_apply = Particle::getB 198 199 //// NN 200 if(particle1->isNucleon() && particle2->is 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; 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 } 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 749 const G4double rChannel=Random::shoot() 750 751 if(elasticCX > rChannel) { 752 isElastic = true; 753 // Elastic PiN channel 754 INCL_DEBUG("PiN interaction: elastic c 755 weight = counterweight; 756 return new PiNElasticChannel(particle1 757 } else if(elasticCX + deltaProductionCX 758 isElastic = false; 759 // PiN -> Delta channel is chosen 760 INCL_DEBUG("PiN interaction: Delta cha 761 weight = counterweight; 762 return new PiNToDeltaChannel(particle1 763 } else if(elasticCX + deltaProductionCX 764 isElastic = false; 765 // PiN -> PiNPi channel is chosen 766 INCL_DEBUG("PiN interaction: one Pion 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; 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 } 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 973 else { 974 INCL_WARN("inconsisten 975 if(twoPiProductionCX>0 976 INCL_WARN("Returni 977 isElastic = false; 978 return new EtaNToP 979 } else if(onePiProduct 980 INCL_WARN("Returni 981 isElastic = false; 982 return new EtaNToP 983 } else { 984 INCL_WARN("Returni 985 isElastic = true; 986 return new EtaNEla 987 } 988 } 989 990 } else if ((particle1->isNucleon() && part 991 //// OmegaN 992 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(); 1369 return NULL; 1370 } 1371 } 1372 1373 void BinaryCollisionAvatar::preInteraction( 1374 isParticle1Spectator = particle1->isTarge 1375 isParticle2Spectator = particle2->isTarge 1376 InteractionAvatar::preInteraction(); 1377 } 1378 1379 void BinaryCollisionAvatar::postInteraction 1380 // Call the postInteraction method of the 1381 // (provides Pauli blocking and enforces 1382 InteractionAvatar::postInteraction(fs); 1383 1384 switch(fs->getValidity()) { 1385 case PauliBlockedFS: 1386 theNucleus->getStore()->getBook().inc 1387 break; 1388 case NoEnergyConservationFS: 1389 case ParticleBelowFermiFS: 1390 case ParticleBelowZeroFS: 1391 break; 1392 case ValidFS: 1393 Book &theBook = theNucleus->getStore( 1394 theBook.incrementAcceptedCollisions() 1395 1396 if(theBook.getAcceptedCollisions() == 1397 // Store time and cross section of 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 } 1419 } 1420 return; 1421 } 1422 1423 std::string BinaryCollisionAvatar::dump() c 1424 std::stringstream ss; 1425 ss << "(avatar " << theTime <<" 'nn-colli 1426 << "(list " << '\n' 1427 << particle1->dump() 1428 << particle2->dump() 1429 << "))" << '\n'; 1430 return ss.str(); 1431 } 1432 1433 } 1434