Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 35 36 #include "globals.hh" 37 38 #include "G4INCLNpiToMissingStrangenessChannel.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 const G4double NpiToMissingStrangenessChannel::angularSlope = 1.; 50 51 NpiToMissingStrangenessChannel::NpiToMissingStrangenessChannel(Particle *p1, Particle *p2) 52 : particle1(p1), particle2(p2) 53 {} 54 55 NpiToMissingStrangenessChannel::~NpiToMissingStrangenessChannel(){} 56 57 void NpiToMissingStrangenessChannel::fillFinalState(FinalState *fs) { 58 59 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2); // MeV /!\/!\/!\. 60 // assert(sqrtS > 2.240); // ! > 2.109 Not supposed to be under 2.244 GeV. 61 62 const G4int iso = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType()); 63 // assert(iso == -3 || iso == -1 || iso == 1 || iso == 3); 64 G4int iso_system = 0; 65 G4int available_iso = 0; 66 G4int nbr_pions = 0; 67 G4int min_pions = 0; 68 G4int max_pions = 0; 69 70 Particle *pion_initial; 71 Particle *nucleon_initial; 72 73 if(particle1->isPion()){ 74 pion_initial = particle1; 75 nucleon_initial = particle2; 76 } 77 else{ 78 pion_initial = particle2; 79 nucleon_initial = particle1; 80 } 81 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion_initial, nucleon_initial); // GeV /!\/!\/!\. 82 83 G4double rdm = Random::shoot(); 84 85 //G4int nbr_particle = 2; 86 87 if(rdm < 0.35){ 88 // Lambda-K chosen 89 nucleon_initial->setType(Lambda); 90 available_iso = 1; 91 min_pions = 3; 92 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Lambda)-ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus)); 93 } 94 else if((iso == 0 && rdm < 0.55) || rdm < 0.5){ 95 // N-K-Kb chosen 96 //nbr_particle++; 97 available_iso = 3; 98 min_pions = 1; 99 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Proton)-2.*ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus)); 100 } 101 else{ 102 // Sigma-K chosen 103 available_iso = 3; 104 min_pions = 3; 105 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(SigmaMinus)-ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus)); 106 } 107 // Gaussian noise + mean value nbr pions fonction energy (choice) 108 G4double intermediaire = min_pions + Random::gauss(2) + std::sqrt(pLab-2.2); 109 nbr_pions = std::min(max_pions,std::max(min_pions,G4int(intermediaire ))); 110 111 available_iso += nbr_pions*2; 112 #ifdef INCLXX_IN_GEANT4_MODE 113 // Erase the parent resonance information of the initial particles 114 particle1->setParentResonancePDGCode(0); 115 particle1->setParentResonanceID(0); 116 particle2->setParentResonancePDGCode(0); 117 particle2->setParentResonanceID(0); 118 #endif 119 //nbr_particle += nbr_pions; 120 121 ParticleList list; 122 ParticleType PionType = PiZero; 123 const ThreeVector &rcol1 = pion_initial->getPosition(); 124 const ThreeVector zero; 125 126 // (pip piz pim) (sp sz sm) (L S Kb) 127 //pip_p 0.63 0.26 0.11 0.73 0.25 0.02 0.42 0.49 0.09 // inital 128 //pip_p 0.54 0.26 0.20 0.73 0.25 0.02 0.42 0.49 0.09 // choice 129 G4bool pip_p = (std::abs(iso) == 3); 130 //piz_p 0.32 0.45 0.23 0.52 0.40 0.08 0.40 0.41 0.19 131 G4bool piz_p = (ParticleTable::getIsospin(pion_initial->getType()) == 0); 132 //pim_p 0.18 0.37 0.45 0.20 0.63 0.17 0.39 0.33 0.28 133 G4bool pim_p = (!pip_p && !piz_p); 134 135 for(Int_t i=0; i<nbr_pions; i++){ 136 Particle *pion = new Particle(PionType,zero,rcol1); 137 if(available_iso-std::abs(iso-iso_system) >= 4){ 138 rdm = Random::shoot(); 139 if((pip_p && rdm < 0.54) || (piz_p && rdm < 0.32) || (pim_p && rdm < 0.45)){ 140 pion->setType(ParticleTable::getPionType(G4int(Math::sign(iso))*2)); //pip/pip/pim 141 iso_system += 2*G4int(Math::sign(iso)); 142 available_iso -= 2; 143 } 144 else if((pip_p && rdm < 0.80) || (piz_p && rdm < 0.77) || (pim_p && rdm < 0.82)){ 145 pion->setType(PiZero); 146 available_iso -= 2; 147 } 148 else{ 149 pion->setType(ParticleTable::getPionType(-G4int(Math::sign(iso))*2)); 150 iso_system -= 2*G4int(Math::sign(iso)); 151 available_iso -= 2; 152 } 153 } 154 else if(available_iso-std::abs(iso-iso_system) == 2){ 155 rdm = Random::shoot(); 156 if((pip_p && rdm < 0.26/0.37 && (Math::sign(iso)*Math::sign(iso-iso_system)+1)) || (pip_p && rdm < 0.26/0.89 && (Math::sign(iso)*Math::sign(iso-iso_system)-1)) || 157 (piz_p && rdm < 0.45/0.68 && (Math::sign(iso)*Math::sign(iso-iso_system)+1)) || (piz_p && rdm < 0.45/0.77 && (Math::sign(iso)*Math::sign(iso-iso_system)-1)) || 158 (pim_p && rdm < 0.37/0.82 && (Math::sign(iso)*Math::sign(iso-iso_system)+1)) || (piz_p && rdm < 0.37/0.55 && (Math::sign(iso)*Math::sign(iso-iso_system)-1))){ 159 pion->setType(PiZero); 160 available_iso -= 2; 161 } 162 else{ 163 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2)); 164 iso_system += Math::sign(iso-iso_system)*2; 165 available_iso -= 2; 166 } 167 } 168 else if(available_iso-std::abs(iso-iso_system) == 0){ 169 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2)); 170 iso_system += Math::sign(iso-iso_system)*2; 171 available_iso -= 2; 172 } 173 else INCL_ERROR("Pion Generation Problem in NpiToMissingStrangenessChannel" << '\n'); 174 list.push_back(pion); 175 } 176 177 if(nucleon_initial->isLambda()){ // K-Lambda 178 // assert(available_iso == 1); 179 pion_initial->setType(ParticleTable::getKaonType(iso-iso_system)); 180 } 181 else if(min_pions == 1){ // N-K-Kb chosen 182 // assert(available_iso == 3); 183 Particle *antikaon = new Particle(KMinus,zero,rcol1); 184 if(std::abs(iso-iso_system) == 3){ 185 pion_initial->setType(ParticleTable::getKaonType((iso-iso_system)/3)); 186 nucleon_initial->setType(ParticleTable::getNucleonType((iso-iso_system)/3)); 187 antikaon->setType(ParticleTable::getAntiKaonType((iso-iso_system)/3)); 188 } 189 else if(std::abs(iso-iso_system) == 1){ // equi-repartition 190 rdm = G4int(Random::shoot()*3.)-1; 191 nucleon_initial->setType(ParticleTable::getNucleonType((G4int(rdm+0.5)*2-1)*(iso_system-iso))); 192 pion_initial->setType(ParticleTable::getKaonType((std::abs(rdm*2)-1)*(iso-iso_system))); 193 antikaon->setType(ParticleTable::getAntiKaonType((G4int(rdm-0.5)*2+1)*(iso-iso_system))); 194 } 195 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n'); 196 list.push_back(antikaon); 197 nbr_pions += 1; // need for addCreatedParticle loop 198 } 199 else{// Sigma-K 200 // assert(available_iso == 3); 201 if(std::abs(iso-iso_system) == 3){ 202 pion_initial->setType(ParticleTable::getKaonType((iso-iso_system)/3)); 203 nucleon_initial->setType(ParticleTable::getSigmaType((iso-iso_system)*2/3)); 204 } 205 else if(std::abs(iso-iso_system) == 1){ 206 rdm = Random::shoot(); 207 if((pip_p && rdm < 0.73) || (piz_p && rdm < 0.32) || (pim_p && rdm < 0.45)){ 208 nucleon_initial->setType(SigmaZero); 209 pion_initial->setType(ParticleTable::getKaonType(iso-iso_system)); 210 } 211 else{ 212 nucleon_initial->setType(ParticleTable::getSigmaType((iso-iso_system)*2)); 213 pion_initial->setType(ParticleTable::getKaonType(iso_system-iso)); 214 } 215 } 216 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n'); 217 } 218 219 list.push_back(pion_initial); 220 list.push_back(nucleon_initial); 221 222 PhaseSpaceGenerator::generateBiased(sqrtS, list, list.size()-1, angularSlope); 223 224 fs->addModifiedParticle(pion_initial); 225 fs->addModifiedParticle(nucleon_initial); 226 for(Int_t i=0; i<nbr_pions; i++) fs->addCreatedParticle(list[i]); 227 228 } 229 } 230