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 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland 31 // Sylvie Leray, CEA-Saclay, France 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 34 #define INCLXX_IN_GEANT4_MODE 1 35 35 36 #include "globals.hh" 36 #include "globals.hh" 37 37 38 #include "G4INCLNNToMissingStrangenessChannel. 38 #include "G4INCLNNToMissingStrangenessChannel.hh" 39 #include "G4INCLKinematicsUtils.hh" 39 #include "G4INCLKinematicsUtils.hh" 40 #include "G4INCLBinaryCollisionAvatar.hh" 40 #include "G4INCLBinaryCollisionAvatar.hh" 41 #include "G4INCLRandom.hh" 41 #include "G4INCLRandom.hh" 42 #include "G4INCLGlobals.hh" 42 #include "G4INCLGlobals.hh" 43 #include "G4INCLLogger.hh" 43 #include "G4INCLLogger.hh" 44 #include <algorithm> 44 #include <algorithm> 45 #include "G4INCLPhaseSpaceGenerator.hh" 45 #include "G4INCLPhaseSpaceGenerator.hh" 46 46 47 namespace G4INCL { 47 namespace G4INCL { 48 48 49 const G4double NNToMissingStrangenessChannel 49 const G4double NNToMissingStrangenessChannel::angularSlope = 1.; 50 50 51 NNToMissingStrangenessChannel::NNToMissingSt 51 NNToMissingStrangenessChannel::NNToMissingStrangenessChannel(Particle *p1, Particle *p2) 52 : particle1(p1), particle2(p2) 52 : particle1(p1), particle2(p2) 53 {} 53 {} 54 54 55 NNToMissingStrangenessChannel::~NNToMissingS 55 NNToMissingStrangenessChannel::~NNToMissingStrangenessChannel(){} 56 56 57 void NNToMissingStrangenessChannel::fillFina 57 void NNToMissingStrangenessChannel::fillFinalState(FinalState *fs) { 58 58 59 const G4double sqrtS = KinematicsUtils::to 59 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2); // MeV 60 const G4double pLab = 0.001*KinematicsUtil 60 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(particle1, particle2); // GeV 61 // assert(sqrtS > 3100); // ! > 3047.34. Not s 61 // assert(sqrtS > 3100); // ! > 3047.34. Not supposed to be under 3,626 GeV. 62 62 63 const G4int iso = ParticleTable::getIsospi 63 const G4int iso = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType()); 64 // assert(iso == -2 || iso == 0 || iso == 2); 64 // assert(iso == -2 || iso == 0 || iso == 2); 65 G4int iso_system = 0; 65 G4int iso_system = 0; 66 G4int available_iso = 0; 66 G4int available_iso = 0; 67 G4int nbr_pions = 0; 67 G4int nbr_pions = 0; 68 G4int min_pions = 0; 68 G4int min_pions = 0; 69 G4int max_pions = 0; 69 G4int max_pions = 0; 70 70 71 G4double rdm = Random::shoot(); 71 G4double rdm = Random::shoot(); 72 72 73 //G4int nbr_particle = 3; << 74 << 75 if(rdm < 0.35){ 73 if(rdm < 0.35){ 76 // N-K-Lambda chosen 74 // N-K-Lambda chosen 77 particle2->setType(Lambda); 75 particle2->setType(Lambda); 78 available_iso = 2; 76 available_iso = 2; 79 min_pions = 3; 77 min_pions = 3; 80 max_pions = G4int((sqrtS-ParticleTable:: 78 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Lambda)-ParticleTable::getINCLMass(Proton)-ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus)); 81 } 79 } 82 else if((iso == 0 && rdm < 0.55) || rdm < 80 else if((iso == 0 && rdm < 0.55) || rdm < 0.5){ 83 // N-N-K-Kb chosen 81 // N-N-K-Kb chosen 84 //nbr_particle++; << 85 available_iso = 4; 82 available_iso = 4; 86 min_pions = 1; 83 min_pions = 1; 87 max_pions = G4int((sqrtS-2.*ParticleTabl 84 max_pions = G4int((sqrtS-2.*ParticleTable::getINCLMass(Proton)-2.*ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus)); 88 } 85 } 89 else{ 86 else{ 90 // N-K-Sigma chosen 87 // N-K-Sigma chosen 91 available_iso = 4; 88 available_iso = 4; 92 min_pions = 3; 89 min_pions = 3; 93 max_pions = G4int((sqrtS-ParticleTable:: 90 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(SigmaMinus)-ParticleTable::getINCLMass(Proton)-ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus)); 94 } 91 } 95 /* Gaussian noise + mean 92 /* Gaussian noise + mean value nbr pions fonction energy (choice)*/ 96 G4double intermediaire = min_pions + Rando 93 G4double intermediaire = min_pions + Random::gauss(2) + std::sqrt(pLab-5); 97 nbr_pions = std::min(max_pions,std::max(mi 94 nbr_pions = std::min(max_pions,std::max(min_pions,G4int(intermediaire ))); 98 95 99 available_iso += nbr_pions*2; 96 available_iso += nbr_pions*2; 100 //nbr_particle += nbr_pions; << 101 97 102 ParticleList list; 98 ParticleList list; 103 ParticleType PionType = PiZero; 99 ParticleType PionType = PiZero; 104 const ThreeVector &rcol1 = particle1->getP 100 const ThreeVector &rcol1 = particle1->getPosition(); 105 const ThreeVector zero; 101 const ThreeVector zero; 106 Particle *kaon = new Particle(KPlus,zero,r 102 Particle *kaon = new Particle(KPlus,zero,rcol1); 107 103 108 for(Int_t i=0; i<nbr_pions; i++){ 104 for(Int_t i=0; i<nbr_pions; i++){ 109 Particle *pion = new Particle(PionType,z 105 Particle *pion = new Particle(PionType,zero,rcol1); 110 if(available_iso-std::abs(iso-iso_system 106 if(available_iso-std::abs(iso-iso_system) >= 4){ // pp(pn) pip:0.40(0.3) piz:0.35(0.4) pim:0.25(0.3) 111 rdm = Random::shoot(); 107 rdm = Random::shoot(); 112 if(((iso == 0) && rdm < 0.3) || ((iso 108 if(((iso == 0) && rdm < 0.3) || ((iso == -2) && rdm < 0.40) || rdm < 0.25){ 113 pion->setType(PiMinus); 109 pion->setType(PiMinus); 114 iso_system -= 2; 110 iso_system -= 2; 115 available_iso -= 2; 111 available_iso -= 2; 116 } 112 } 117 else if(((iso == 0) && rdm < 0.7) || ( 113 else if(((iso == 0) && rdm < 0.7) || ((iso == -2) && rdm < 0.75) || rdm < 0.60){ 118 pion->setType(PiZero); 114 pion->setType(PiZero); 119 available_iso -= 2; 115 available_iso -= 2; 120 } 116 } 121 else{ 117 else{ 122 pion->setType(PiPlus); 118 pion->setType(PiPlus); 123 iso_system += 2; 119 iso_system += 2; 124 available_iso -= 2; 120 available_iso -= 2; 125 } 121 } 126 } 122 } 127 else if(available_iso-std::abs(iso-iso_s 123 else if(available_iso-std::abs(iso-iso_system) == 2){ 128 rdm = Random::shoot(); 124 rdm = Random::shoot(); 129 // pn pp too hig 125 // pn pp too high (nn too low) -> PiMinus(PiPlus) pp too low (nn too high) -> PiPlus(PiMinus) 130 if(((iso == 0) && (rdm*0.7 < 0.3)) || 126 if(((iso == 0) && (rdm*0.7 < 0.3)) || ((rdm*0.60 < 0.25) && (Math::sign(iso-iso_system)*2-iso != 0)) || ((rdm*0.75 < 0.40) && (Math::sign(iso-iso_system)*2+iso != 0) && (iso != 0))){ 131 pion->setType(ParticleTable::getPion 127 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2)); 132 iso_system += Math::sign(iso-iso_sys 128 iso_system += Math::sign(iso-iso_system)*2; 133 available_iso -= 2; 129 available_iso -= 2; 134 } 130 } 135 else{ 131 else{ 136 PionType = PiZero; 132 PionType = PiZero; 137 available_iso -= 2; 133 available_iso -= 2; 138 } 134 } 139 } 135 } 140 else if(available_iso-std::abs(iso-iso_s 136 else if(available_iso-std::abs(iso-iso_system) == 0){ 141 pion->setType(ParticleTable::getPionTy 137 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2)); 142 iso_system += Math::sign(iso-iso_syste 138 iso_system += Math::sign(iso-iso_system)*2; 143 available_iso -= 2; 139 available_iso -= 2; 144 } 140 } 145 else INCL_ERROR("Pion Generation Problem 141 else INCL_ERROR("Pion Generation Problem in NNToMissingStrangenessChannel" << '\n'); 146 list.push_back(pion); 142 list.push_back(pion); 147 } 143 } 148 144 149 if(particle2->isLambda()){ // N-K-Lambda 145 if(particle2->isLambda()){ // N-K-Lambda 150 // assert(available_iso == 2); 146 // assert(available_iso == 2); 151 if(std::abs(iso-iso_system) == 2){ 147 if(std::abs(iso-iso_system) == 2){ 152 particle1->setType(ParticleTable::getN 148 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2)); 153 kaon->setType(ParticleTable::getKaonTy 149 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/2)); 154 } 150 } 155 else if(std::abs(iso-iso_system) == 0){ 151 else if(std::abs(iso-iso_system) == 0){ 156 rdm = G4int(Random::shoot()*2.)*2-1; 152 rdm = G4int(Random::shoot()*2.)*2-1; 157 particle1->setType(ParticleTable::getN 153 particle1->setType(ParticleTable::getNucleonType(G4int(rdm))); 158 kaon->setType(ParticleTable::getKaonTy 154 kaon->setType(ParticleTable::getKaonType(G4int(-rdm))); 159 } 155 } 160 else INCL_ERROR("Isospin non-conservatio 156 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n'); 161 } 157 } 162 else if(min_pions == 1){ // N-N-K-Kb chose 158 else if(min_pions == 1){ // N-N-K-Kb chosen 163 // assert(available_iso == 4); 159 // assert(available_iso == 4); 164 Particle *antikaon = new Particle(KMinus 160 Particle *antikaon = new Particle(KMinus,zero,rcol1); 165 if(std::abs(iso-iso_system) == 4){ 161 if(std::abs(iso-iso_system) == 4){ 166 particle1->setType(ParticleTable::getN 162 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/4)); 167 particle2->setType(ParticleTable::getN 163 particle2->setType(ParticleTable::getNucleonType((iso-iso_system)/4)); 168 kaon->setType(ParticleTable::getKaonTy 164 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/4)); 169 antikaon->setType(ParticleTable::getAn 165 antikaon->setType(ParticleTable::getAntiKaonType((iso-iso_system)/4)); 170 } 166 } 171 else if(std::abs(iso-iso_system) == 2){ 167 else if(std::abs(iso-iso_system) == 2){ // choice: kaon type free, nucleon type fixed 172 rdm = G4int(Random::shoot()*2.)*2-1; 168 rdm = G4int(Random::shoot()*2.)*2-1; 173 particle1->setType(ParticleTable::getN 169 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2)); 174 particle2->setType(ParticleTable::getN 170 particle2->setType(ParticleTable::getNucleonType((iso-iso_system)/2)); 175 kaon->setType(ParticleTable::getKaonTy 171 kaon->setType(ParticleTable::getKaonType(G4int(rdm))); 176 antikaon->setType(ParticleTable::getAn 172 antikaon->setType(ParticleTable::getAntiKaonType(G4int(-rdm))); 177 } 173 } 178 else if(std::abs(iso-iso_system) == 0){ 174 else if(std::abs(iso-iso_system) == 0){ // particle1 3/4 proton, 1/4 neutron; particle2 1/4 proton, 3/4 neutron 179 rdm = G4int(Random::shoot()*2.)*2-1; 175 rdm = G4int(Random::shoot()*2.)*2-1; 180 G4double rdm2 = G4int(Random::shoot()* 176 G4double rdm2 = G4int(Random::shoot()*2.)*2-1; 181 particle1->setType(ParticleTable::getN 177 particle1->setType(ParticleTable::getNucleonType(std::max(rdm,rdm2))); 182 particle2->setType(ParticleTable::getN 178 particle2->setType(ParticleTable::getNucleonType(std::min(rdm,rdm2))); 183 kaon->setType(ParticleTable::getKaonTy 179 kaon->setType(ParticleTable::getKaonType(-G4int(rdm))); 184 antikaon->setType(ParticleTable::getAn 180 antikaon->setType(ParticleTable::getAntiKaonType(-G4int(rdm2))); 185 } 181 } 186 else INCL_ERROR("Isospin non-conservatio 182 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n'); 187 list.push_back(antikaon); 183 list.push_back(antikaon); 188 nbr_pions += 1; // need for addCreatedPa 184 nbr_pions += 1; // need for addCreatedParticle loop 189 } 185 } 190 else{// N-K-Sigma 186 else{// N-K-Sigma 191 // assert(available_iso == 4); 187 // assert(available_iso == 4); 192 if(std::abs(iso-iso_system) == 4){ 188 if(std::abs(iso-iso_system) == 4){ 193 particle1->setType(ParticleTable::getN 189 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/4)); 194 particle2->setType(ParticleTable::getS 190 particle2->setType(ParticleTable::getSigmaType((iso-iso_system)/2)); 195 kaon->setType(ParticleTable::getKaonTy 191 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/4)); 196 } 192 } 197 else if(std::abs(iso-iso_system) == 2){ 193 else if(std::abs(iso-iso_system) == 2){ // choice: sigma quasi-free, kaon type semi-free, nucleon type fixed 198 rdm = Random::shoot(); 194 rdm = Random::shoot(); 199 // pp(pn) sp:0.42(0.23) sz:0.51(0.54) 195 // pp(pn) sp:0.42(0.23) sz:0.51(0.54) sm:0.07(0.23) 200 if(((iso == 0) && (rdm*0.77 < 0.23)) | 196 if(((iso == 0) && (rdm*0.77 < 0.23)) || ((rdm*0.58 < 0.07) && (Math::sign(iso-iso_system)*2-iso != 0)) || ((rdm*0.93 < 0.42) && (Math::sign(iso-iso_system)*2+iso != 0) && (iso != 0))){ 201 particle2->setType(ParticleTable::ge 197 particle2->setType(ParticleTable::getSigmaType(Math::sign(iso-iso_system)*2)); 202 rdm = G4int(Random::shoot()*2.)*2-1; 198 rdm = G4int(Random::shoot()*2.)*2-1; 203 particle1->setType(ParticleTable::ge 199 particle1->setType(ParticleTable::getNucleonType(G4int(rdm))); 204 kaon->setType(ParticleTable::getKaon 200 kaon->setType(ParticleTable::getKaonType(-G4int(rdm))); 205 } 201 } 206 else{ 202 else{ 207 particle2->setType(SigmaZero); 203 particle2->setType(SigmaZero); 208 particle1->setType(ParticleTable::ge 204 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2)); 209 kaon->setType(ParticleTable::getKaon 205 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/2)); 210 } 206 } 211 } 207 } 212 else if(std::abs(iso-iso_system) == 0){ 208 else if(std::abs(iso-iso_system) == 0){ // choice: sigma free, kaontype semi-free, nucleon type fixed 213 if(((iso == 0) && rdm < 0.23) || ((iso 209 if(((iso == 0) && rdm < 0.23) || ((iso == 2) && rdm < 0.42) || rdm < 0.07){ 214 particle2->setType(SigmaPlus); 210 particle2->setType(SigmaPlus); 215 particle1->setType(Neutron); 211 particle1->setType(Neutron); 216 kaon->setType(KZero); 212 kaon->setType(KZero); 217 } 213 } 218 else if(((iso == 0) && rdm < 0.77) || 214 else if(((iso == 0) && rdm < 0.77) || ((iso == 2) && rdm < 0.93) || rdm < 0.58){ 219 particle2->setType(SigmaZero); 215 particle2->setType(SigmaZero); 220 rdm = G4int(Random::shoot()*2.)*2-1; 216 rdm = G4int(Random::shoot()*2.)*2-1; 221 particle1->setType(ParticleTable::ge 217 particle1->setType(ParticleTable::getNucleonType(G4int(rdm))); 222 kaon->setType(ParticleTable::getKaon 218 kaon->setType(ParticleTable::getKaonType(-G4int(rdm))); 223 } 219 } 224 else{ 220 else{ 225 particle2->setType(SigmaMinus); 221 particle2->setType(SigmaMinus); 226 particle1->setType(Proton); 222 particle1->setType(Proton); 227 kaon->setType(KPlus); 223 kaon->setType(KPlus); 228 } 224 } 229 } 225 } 230 else INCL_ERROR("Isospin non-conservatio 226 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n'); 231 } 227 } 232 228 233 list.push_back(particle1); 229 list.push_back(particle1); 234 list.push_back(particle2); 230 list.push_back(particle2); 235 list.push_back(kaon); 231 list.push_back(kaon); 236 232 237 PhaseSpaceGenerator::generateBiased(sqrtS, 233 PhaseSpaceGenerator::generateBiased(sqrtS, list, list.size()-3, angularSlope); 238 234 239 fs->addModifiedParticle(particle1); 235 fs->addModifiedParticle(particle1); 240 fs->addModifiedParticle(particle2); 236 fs->addModifiedParticle(particle2); 241 fs->addCreatedParticle(kaon); 237 fs->addCreatedParticle(kaon); 242 for(Int_t i=0; i<nbr_pions; i++) fs->addCr 238 for(Int_t i=0; i<nbr_pions; i++) fs->addCreatedParticle(list[i]); 243 239 244 } 240 } 245 } 241 } 246 242