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 #include "G4INCLNNbarToLLbarChannel.hh" 39 #include "G4INCLKinematicsUtils.hh" 40 #include "G4INCLBinaryCollisionAvatar.hh" 41 #include "G4INCLRandom.hh" 42 #include "G4INCLGlobals.hh" 43 #include "G4INCLLogger.hh" 44 #include <algorithm> 45 #include "G4INCLPhaseSpaceGenerator.hh" 46 47 namespace G4INCL { 48 49 NNbarToLLbarChannel::NNbarToLLbarChannel(Par 50 : particle1(p1), particle2(p2) 51 {} 52 53 NNbarToLLbarChannel::~NNbarToLLbarChannel(){ 54 55 void NNbarToLLbarChannel::fillFinalState(Fin 56 // this channel include all states with la 57 58 //brief ppbar 59 // p pbar -> l lbar (BFMM 121) 60 // ppbar -> l lbar pi0 (BFMM 113) 61 // ppbar -> splus pim lbar || sminusba 62 // ppbar -> sminus pip lbar || splusba 63 // ppbar -> sp spbar (BFMM 139) 64 // ppbar -> sm smbar (BFMM 149) 65 // ppbar -> szero szerobar (BFMM 144) 66 // ppbar -> ximinus ximinusbar (BFMM 1 67 // ppbar -> szero lbar || szerobar l ( 68 // 69 // 70 //brief npbar 71 // n pbar -> l lbar pi- (BFMM 487) 72 // n pbar -> l sbarplus || lbar sminus 73 // 74 // 75 //brief nnbar 76 // all same as for ppbar 77 // 78 // 79 //brief pnbar 80 // p nbar -> l lbar pi+ (same as BFMM 81 // p nbar -> l sbarminus || lbar splus 82 // 83 84 Particle *nucleon; 85 Particle *antinucleon; 86 87 if(particle1->isNucleon()){ 88 nucleon = particle1; 89 antinucleon = particle2; 90 } 91 else{ 92 nucleon = particle2; 93 antinucleon = particle1; 94 } 95 96 const G4double plab = 0.001*KinematicsUtil 97 // ppbar cross sections 98 99 const std::vector<G4double> BFMM121 = {2.3 100 //const G4double Eth_PPbar_LLbar = 1.437; 101 const std::vector<G4double> BFMM113 = {-0. 102 //const G4double Eth_PPbar_LLbar_pi0 = 1.8 103 const std::vector<G4double> BFMM139 = {0.1 104 //const G4double Eth_PPbar_SpSpbar = 1.851 105 const std::vector<G4double> BFMM149 = {1.8 106 //const G4double Eth_PPbar_SmSmbar = 1.896 107 const std::vector<G4double> BFMM136 = {1.7 108 //const G4double Eth_PPbar_SpLbar_pim = 2. 109 const std::vector<G4double> BFMM146 = {1.0 110 //const G4double Eth_PPbar_SmLbar_pip = 2. 111 const std::vector<G4double> BFMM143 = {0.6 112 113 114 115 //const G4double Eth_PPbar_Szero_Lbar = 1. 116 //fixed due to limited data 117 G4double BFMM144; 118 if(plab > 2.0) BFMM144 = 0.008; //sigmazer 119 else BFMM144 = 0.0; 120 G4double BFMM101; 121 if(plab > 2.8) BFMM101 = 0.002; //ximinus 122 else BFMM101 = 0.0; 123 124 // npbar cross sections (fixed due to limi 125 G4double BFMM487; 126 if(plab > 2.1) BFMM487 = 0.048; //llbar pi 127 else BFMM487 = 0.0; 128 G4double BFMM488; 129 if(plab > 2.0) BFMM488 = 0.139; //lsigmami 130 else BFMM488 = 0.0; 131 132 const G4double sqrtS = KinematicsUtils::to 133 const G4double totalppbar = KinematicsUtil 134 +KinematicsUtils::compute_xs(BFMM139, plab 135 +KinematicsUtils::compute_xs(BFMM146, plab 136 +KinematicsUtils::compute_xs(BFMM121, plab 137 +BFMM144 +BFMM101; 138 const G4double totalpnbar = BFMM487 + BFMM 139 const G4double rdm = Random::shoot(); 140 141 G4bool thirdparticle = false; //set true i 142 ParticleType PionType; 143 //setting types of new particles 144 if(nucleon->getType()==Proton){ 145 if(antinucleon->getType()==antiProton){ 146 if(rdm*totalppbar < KinematicsUtils::c 147 nucleon->setType(Lambda); 148 antinucleon->setType(antiLambda); 149 } 150 else if(rdm*totalppbar < KinematicsUti 151 nucleon->setType(SigmaZero); 152 antinucleon->setType(antiSigmaZero); 153 } 154 else if(rdm*totalppbar < KinematicsUti 155 nucleon->setType(XiMinus); 156 antinucleon->setType(antiXiMinus); 157 } 158 else if(rdm*totalppbar < KinematicsUti 159 +KinematicsUtils::compute_xs(BFMM113, 160 nucleon->setType(Lambda); 161 antinucleon->setType(antiLambda); 162 thirdparticle = true; 163 PionType = PiZero; 164 } 165 else if(rdm*totalppbar < KinematicsUti 166 +KinematicsUtils::compute_xs(BFMM113, 167 G4double rdm2 = Random::shoot(); 168 if(rdm2 > 0.5){ 169 nucleon->setType(SigmaPlus); 170 antinucleon->setType(antiLambda); 171 thirdparticle = true; 172 PionType = PiMinus; 173 } 174 else{ 175 nucleon->setType(antiSigmaMinus); 176 antinucleon->setType(Lambda); 177 thirdparticle = true; 178 PionType = PiMinus; 179 } 180 } 181 else if(rdm*totalppbar < KinematicsUti 182 +KinematicsUtils::compute_xs(BFMM113, 183 +KinematicsUtils::compute_xs(BFMM146, 184 G4double rdm2 = Random::shoot(); 185 if(rdm2 > 0.5){ 186 nucleon->setType(SigmaMinus); 187 antinucleon->setType(antiLambda); 188 thirdparticle = true; 189 PionType = PiPlus; 190 } 191 else{ 192 nucleon->setType(antiSigmaPlus); 193 antinucleon->setType(Lambda); 194 thirdparticle = true; 195 PionType = PiPlus; 196 } 197 } 198 else if(rdm*totalppbar < KinematicsUti 199 +KinematicsUtils::compute_xs(BFMM113, 200 +KinematicsUtils::compute_xs(BFMM146, 201 G4double rdm2 = Random::shoot(); 202 if(rdm2 > 0.5){ 203 nucleon->setType(SigmaZero); 204 antinucleon->setType(antiLambda); 205 } 206 else{ 207 nucleon->setType(antiSigmaZero); 208 antinucleon->setType(Lambda); 209 } 210 } 211 else if(rdm*totalppbar < KinematicsUti 212 +KinematicsUtils::compute_xs(BFMM113, 213 +KinematicsUtils::compute_xs(BFMM146, 214 +KinematicsUtils::compute_xs(BFMM139, 215 nucleon->setType(SigmaPlus); 216 antinucleon->setType(antiSigmaPlus); 217 } 218 else if(rdm*totalppbar < KinematicsUti 219 +KinematicsUtils::compute_xs(std::move 220 +KinematicsUtils::compute_xs(std::move 221 +KinematicsUtils::compute_xs(std::move 222 nucleon->setType(SigmaMinus); 223 antinucleon->setType(antiSigmaMinus) 224 } 225 else{ 226 INCL_ERROR("out of total ppbar sum i 227 } 228 } 229 else{ //pnbar case charge +1 230 if(rdm*totalpnbar < BFMM488){ 231 G4double rdm2 = Random::shoot(); 232 if(rdm2 > 0.5){ 233 nucleon->setType(Lambda); 234 antinucleon->setType(antiSigmaMinu 235 } 236 else{ 237 nucleon->setType(antiLambda); 238 antinucleon->setType(SigmaPlus); / 239 } 240 } 241 else{ 242 nucleon->setType(Lambda); 243 antinucleon->setType(antiLambda); 244 thirdparticle = true; 245 PionType = PiPlus; 246 } 247 } 248 } 249 else{ // neutron 250 if(antinucleon->getType()==antiNeutron){ 251 if(rdm*totalppbar < KinematicsUtils::c 252 nucleon->setType(Lambda); 253 antinucleon->setType(antiLambda); 254 } 255 else if(rdm*totalppbar < KinematicsUti 256 nucleon->setType(SigmaZero); 257 antinucleon->setType(antiSigmaZero); 258 } 259 else if(rdm*totalppbar < KinematicsUti 260 nucleon->setType(XiMinus); 261 antinucleon->setType(antiXiMinus); 262 } 263 else if(rdm*totalppbar < KinematicsUti 264 +KinematicsUtils::compute_xs(BFMM113, 265 nucleon->setType(Lambda); 266 antinucleon->setType(antiLambda); 267 thirdparticle = true; 268 PionType = PiZero; 269 } 270 else if(rdm*totalppbar < KinematicsUti 271 +KinematicsUtils::compute_xs(BFMM113, 272 G4double rdm2 = Random::shoot(); 273 if(rdm2 > 0.5){ 274 nucleon->setType(SigmaPlus); 275 antinucleon->setType(antiLambda); 276 thirdparticle = true; 277 PionType = PiMinus; 278 } 279 else{ 280 nucleon->setType(antiSigmaMinus); 281 antinucleon->setType(Lambda); 282 thirdparticle = true; 283 PionType = PiMinus; 284 } 285 } 286 else if(rdm*totalppbar < KinematicsUti 287 +KinematicsUtils::compute_xs(BFMM113, 288 +KinematicsUtils::compute_xs(BFMM146, 289 G4double rdm2 = Random::shoot(); 290 if(rdm2 > 0.5){ 291 nucleon->setType(SigmaMinus); //ch 292 antinucleon->setType(antiLambda); 293 thirdparticle = true; 294 PionType = PiPlus; 295 } 296 else{ 297 nucleon->setType(antiSigmaPlus); / 298 antinucleon->setType(Lambda); 299 thirdparticle = true; 300 PionType = PiPlus; 301 } 302 } 303 else if(rdm*totalppbar < KinematicsUti 304 +KinematicsUtils::compute_xs(BFMM113, 305 +KinematicsUtils::compute_xs(BFMM146, 306 G4double rdm2 = Random::shoot(); 307 if(rdm2 > 0.5){ 308 nucleon->setType(SigmaZero); 309 antinucleon->setType(antiLambda); 310 } 311 else{ 312 nucleon->setType(antiSigmaZero); 313 antinucleon->setType(Lambda); 314 } 315 } 316 else if(rdm*totalppbar < KinematicsUti 317 +KinematicsUtils::compute_xs(BFMM113, 318 +KinematicsUtils::compute_xs(BFMM146, 319 +KinematicsUtils::compute_xs(BFMM139, 320 nucleon->setType(SigmaPlus); 321 antinucleon->setType(antiSigmaPlus); 322 } 323 else if(rdm*totalppbar < KinematicsUti 324 +KinematicsUtils::compute_xs(std::move 325 +KinematicsUtils::compute_xs(std::move 326 +KinematicsUtils::compute_xs(std::move 327 nucleon->setType(SigmaMinus); 328 antinucleon->setType(antiSigmaMinus) 329 } 330 else{ 331 INCL_ERROR("out of total nnbar sum i 332 } 333 } 334 else{ //npbar case charge -1 335 if(rdm*totalpnbar < BFMM488){ 336 G4double rdm2 = Random::shoot(); 337 if(rdm2 > 0.5){ 338 nucleon->setType(Lambda); 339 antinucleon->setType(antiSigmaPlus 340 } 341 else{ 342 nucleon->setType(antiLambda); 343 antinucleon->setType(SigmaMinus); 344 } 345 } 346 else{ 347 nucleon->setType(Lambda); 348 antinucleon->setType(antiLambda); 349 thirdparticle = true; 350 PionType = PiMinus; 351 } 352 } 353 } 354 355 //now assigning momentum to the final part 356 357 if(thirdparticle){ //three particles 358 ParticleList list; 359 list.push_back(nucleon); 360 list.push_back(antinucleon); 361 const ThreeVector &rcol = nucleon->getPo 362 const ThreeVector zero; 363 Particle *pion = new Particle(PionType,z 364 list.push_back(pion); 365 366 PhaseSpaceGenerator::generate(sqrtS, lis 367 368 fs->addModifiedParticle(nucleon); 369 fs->addModifiedParticle(antinucleon); 370 fs->addCreatedParticle(pion); 371 } 372 else{//only two particles 373 G4double mn=nucleon->getMass(); 374 G4double my=antinucleon->getMass(); 375 376 G4double ey=(sqrtS*sqrtS+my*my-mn*mn)/(2 377 G4double en=std::sqrt(ey*ey-my*my+mn*mn) 378 nucleon->setEnergy(en); 379 antinucleon->setEnergy(ey); 380 G4double py=std::sqrt(ey*ey-my*my); 381 382 ThreeVector mom_antinucleon = Random::no 383 384 antinucleon->setMomentum(mom_antinucleon 385 nucleon->setMomentum(-mom_antinucleon); 386 387 fs->addModifiedParticle(nucleon); 388 fs->addModifiedParticle(antinucleon); 389 } 390 391 } 392 } 393