| 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; << 73 G4int nbr_particle = 3;
74 74
75 if(rdm < 0.35){ 75 if(rdm < 0.35){
76 // N-K-Lambda chosen 76 // N-K-Lambda chosen
77 particle2->setType(Lambda); 77 particle2->setType(Lambda);
78 available_iso = 2; 78 available_iso = 2;
79 min_pions = 3; 79 min_pions = 3;
80 max_pions = G4int((sqrtS-ParticleTable:: 80 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Lambda)-ParticleTable::getINCLMass(Proton)-ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus));
81 } 81 }
82 else if((iso == 0 && rdm < 0.55) || rdm < 82 else if((iso == 0 && rdm < 0.55) || rdm < 0.5){
83 // N-N-K-Kb chosen 83 // N-N-K-Kb chosen
84 //nbr_particle++; << 84 nbr_particle++;
85 available_iso = 4; 85 available_iso = 4;
86 min_pions = 1; 86 min_pions = 1;
87 max_pions = G4int((sqrtS-2.*ParticleTabl 87 max_pions = G4int((sqrtS-2.*ParticleTable::getINCLMass(Proton)-2.*ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus));
88 } 88 }
89 else{ 89 else{
90 // N-K-Sigma chosen 90 // N-K-Sigma chosen
91 available_iso = 4; 91 available_iso = 4;
92 min_pions = 3; 92 min_pions = 3;
93 max_pions = G4int((sqrtS-ParticleTable:: 93 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(SigmaMinus)-ParticleTable::getINCLMass(Proton)-ParticleTable::getINCLMass(KZero))/ParticleTable::getINCLMass(PiPlus));
94 } 94 }
95 /* Gaussian noise + mean 95 /* Gaussian noise + mean value nbr pions fonction energy (choice)*/
96 G4double intermediaire = min_pions + Rando 96 G4double intermediaire = min_pions + Random::gauss(2) + std::sqrt(pLab-5);
97 nbr_pions = std::min(max_pions,std::max(mi 97 nbr_pions = std::min(max_pions,std::max(min_pions,G4int(intermediaire )));
98 98
99 available_iso += nbr_pions*2; 99 available_iso += nbr_pions*2;
100 //nbr_particle += nbr_pions; << 100 nbr_particle += nbr_pions;
101 101
102 ParticleList list; 102 ParticleList list;
103 ParticleType PionType = PiZero; 103 ParticleType PionType = PiZero;
104 const ThreeVector &rcol1 = particle1->getP 104 const ThreeVector &rcol1 = particle1->getPosition();
105 const ThreeVector zero; 105 const ThreeVector zero;
106 Particle *kaon = new Particle(KPlus,zero,r 106 Particle *kaon = new Particle(KPlus,zero,rcol1);
107 107
108 for(Int_t i=0; i<nbr_pions; i++){ 108 for(Int_t i=0; i<nbr_pions; i++){
109 Particle *pion = new Particle(PionType,z 109 Particle *pion = new Particle(PionType,zero,rcol1);
110 if(available_iso-std::abs(iso-iso_system 110 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(); 111 rdm = Random::shoot();
112 if(((iso == 0) && rdm < 0.3) || ((iso 112 if(((iso == 0) && rdm < 0.3) || ((iso == -2) && rdm < 0.40) || rdm < 0.25){
113 pion->setType(PiMinus); 113 pion->setType(PiMinus);
114 iso_system -= 2; 114 iso_system -= 2;
115 available_iso -= 2; 115 available_iso -= 2;
116 } 116 }
117 else if(((iso == 0) && rdm < 0.7) || ( 117 else if(((iso == 0) && rdm < 0.7) || ((iso == -2) && rdm < 0.75) || rdm < 0.60){
118 pion->setType(PiZero); 118 pion->setType(PiZero);
119 available_iso -= 2; 119 available_iso -= 2;
120 } 120 }
121 else{ 121 else{
122 pion->setType(PiPlus); 122 pion->setType(PiPlus);
123 iso_system += 2; 123 iso_system += 2;
124 available_iso -= 2; 124 available_iso -= 2;
125 } 125 }
126 } 126 }
127 else if(available_iso-std::abs(iso-iso_s 127 else if(available_iso-std::abs(iso-iso_system) == 2){
128 rdm = Random::shoot(); 128 rdm = Random::shoot();
129 // pn pp too hig << 129 // 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)) || 130 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 131 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
132 iso_system += Math::sign(iso-iso_sys 132 iso_system += Math::sign(iso-iso_system)*2;
133 available_iso -= 2; 133 available_iso -= 2;
134 } 134 }
135 else{ 135 else{
136 PionType = PiZero; 136 PionType = PiZero;
137 available_iso -= 2; 137 available_iso -= 2;
138 } 138 }
139 } 139 }
140 else if(available_iso-std::abs(iso-iso_s 140 else if(available_iso-std::abs(iso-iso_system) == 0){
141 pion->setType(ParticleTable::getPionTy 141 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
142 iso_system += Math::sign(iso-iso_syste 142 iso_system += Math::sign(iso-iso_system)*2;
143 available_iso -= 2; 143 available_iso -= 2;
144 } 144 }
145 else INCL_ERROR("Pion Generation Problem 145 else INCL_ERROR("Pion Generation Problem in NNToMissingStrangenessChannel" << '\n');
146 list.push_back(pion); 146 list.push_back(pion);
147 } 147 }
148 148
149 if(particle2->isLambda()){ // N-K-Lambda 149 if(particle2->isLambda()){ // N-K-Lambda
150 // assert(available_iso == 2); 150 // assert(available_iso == 2);
151 if(std::abs(iso-iso_system) == 2){ 151 if(std::abs(iso-iso_system) == 2){
152 particle1->setType(ParticleTable::getN 152 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
153 kaon->setType(ParticleTable::getKaonTy 153 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/2));
154 } 154 }
155 else if(std::abs(iso-iso_system) == 0){ 155 else if(std::abs(iso-iso_system) == 0){
156 rdm = G4int(Random::shoot()*2.)*2-1; 156 rdm = G4int(Random::shoot()*2.)*2-1;
157 particle1->setType(ParticleTable::getN 157 particle1->setType(ParticleTable::getNucleonType(G4int(rdm)));
158 kaon->setType(ParticleTable::getKaonTy 158 kaon->setType(ParticleTable::getKaonType(G4int(-rdm)));
159 } 159 }
160 else INCL_ERROR("Isospin non-conservatio 160 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
161 } 161 }
162 else if(min_pions == 1){ // N-N-K-Kb chose 162 else if(min_pions == 1){ // N-N-K-Kb chosen
163 // assert(available_iso == 4); 163 // assert(available_iso == 4);
164 Particle *antikaon = new Particle(KMinus 164 Particle *antikaon = new Particle(KMinus,zero,rcol1);
165 if(std::abs(iso-iso_system) == 4){ 165 if(std::abs(iso-iso_system) == 4){
166 particle1->setType(ParticleTable::getN 166 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/4));
167 particle2->setType(ParticleTable::getN 167 particle2->setType(ParticleTable::getNucleonType((iso-iso_system)/4));
168 kaon->setType(ParticleTable::getKaonTy 168 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/4));
169 antikaon->setType(ParticleTable::getAn 169 antikaon->setType(ParticleTable::getAntiKaonType((iso-iso_system)/4));
170 } 170 }
171 else if(std::abs(iso-iso_system) == 2){ 171 else if(std::abs(iso-iso_system) == 2){ // choice: kaon type free, nucleon type fixed
172 rdm = G4int(Random::shoot()*2.)*2-1; 172 rdm = G4int(Random::shoot()*2.)*2-1;
173 particle1->setType(ParticleTable::getN 173 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
174 particle2->setType(ParticleTable::getN 174 particle2->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
175 kaon->setType(ParticleTable::getKaonTy 175 kaon->setType(ParticleTable::getKaonType(G4int(rdm)));
176 antikaon->setType(ParticleTable::getAn 176 antikaon->setType(ParticleTable::getAntiKaonType(G4int(-rdm)));
177 } 177 }
178 else if(std::abs(iso-iso_system) == 0){ 178 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; 179 rdm = G4int(Random::shoot()*2.)*2-1;
180 G4double rdm2 = G4int(Random::shoot()* 180 G4double rdm2 = G4int(Random::shoot()*2.)*2-1;
181 particle1->setType(ParticleTable::getN 181 particle1->setType(ParticleTable::getNucleonType(std::max(rdm,rdm2)));
182 particle2->setType(ParticleTable::getN 182 particle2->setType(ParticleTable::getNucleonType(std::min(rdm,rdm2)));
183 kaon->setType(ParticleTable::getKaonTy 183 kaon->setType(ParticleTable::getKaonType(-G4int(rdm)));
184 antikaon->setType(ParticleTable::getAn 184 antikaon->setType(ParticleTable::getAntiKaonType(-G4int(rdm2)));
185 } 185 }
186 else INCL_ERROR("Isospin non-conservatio 186 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
187 list.push_back(antikaon); 187 list.push_back(antikaon);
188 nbr_pions += 1; // need for addCreatedPa 188 nbr_pions += 1; // need for addCreatedParticle loop
189 } 189 }
190 else{// N-K-Sigma 190 else{// N-K-Sigma
191 // assert(available_iso == 4); 191 // assert(available_iso == 4);
192 if(std::abs(iso-iso_system) == 4){ 192 if(std::abs(iso-iso_system) == 4){
193 particle1->setType(ParticleTable::getN 193 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/4));
194 particle2->setType(ParticleTable::getS 194 particle2->setType(ParticleTable::getSigmaType((iso-iso_system)/2));
195 kaon->setType(ParticleTable::getKaonTy 195 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/4));
196 } 196 }
197 else if(std::abs(iso-iso_system) == 2){ 197 else if(std::abs(iso-iso_system) == 2){ // choice: sigma quasi-free, kaon type semi-free, nucleon type fixed
198 rdm = Random::shoot(); 198 rdm = Random::shoot();
199 // pp(pn) sp:0.42(0.23) sz:0.51(0.54) 199 // 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)) | 200 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 201 particle2->setType(ParticleTable::getSigmaType(Math::sign(iso-iso_system)*2));
202 rdm = G4int(Random::shoot()*2.)*2-1; 202 rdm = G4int(Random::shoot()*2.)*2-1;
203 particle1->setType(ParticleTable::ge 203 particle1->setType(ParticleTable::getNucleonType(G4int(rdm)));
204 kaon->setType(ParticleTable::getKaon 204 kaon->setType(ParticleTable::getKaonType(-G4int(rdm)));
205 } 205 }
206 else{ 206 else{
207 particle2->setType(SigmaZero); 207 particle2->setType(SigmaZero);
208 particle1->setType(ParticleTable::ge 208 particle1->setType(ParticleTable::getNucleonType((iso-iso_system)/2));
209 kaon->setType(ParticleTable::getKaon 209 kaon->setType(ParticleTable::getKaonType((iso-iso_system)/2));
210 } 210 }
211 } 211 }
212 else if(std::abs(iso-iso_system) == 0){ 212 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 213 if(((iso == 0) && rdm < 0.23) || ((iso == 2) && rdm < 0.42) || rdm < 0.07){
214 particle2->setType(SigmaPlus); 214 particle2->setType(SigmaPlus);
215 particle1->setType(Neutron); 215 particle1->setType(Neutron);
216 kaon->setType(KZero); 216 kaon->setType(KZero);
217 } 217 }
218 else if(((iso == 0) && rdm < 0.77) || 218 else if(((iso == 0) && rdm < 0.77) || ((iso == 2) && rdm < 0.93) || rdm < 0.58){
219 particle2->setType(SigmaZero); 219 particle2->setType(SigmaZero);
220 rdm = G4int(Random::shoot()*2.)*2-1; 220 rdm = G4int(Random::shoot()*2.)*2-1;
221 particle1->setType(ParticleTable::ge 221 particle1->setType(ParticleTable::getNucleonType(G4int(rdm)));
222 kaon->setType(ParticleTable::getKaon 222 kaon->setType(ParticleTable::getKaonType(-G4int(rdm)));
223 } 223 }
224 else{ 224 else{
225 particle2->setType(SigmaMinus); 225 particle2->setType(SigmaMinus);
226 particle1->setType(Proton); 226 particle1->setType(Proton);
227 kaon->setType(KPlus); 227 kaon->setType(KPlus);
228 } 228 }
229 } 229 }
230 else INCL_ERROR("Isospin non-conservatio 230 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
231 } 231 }
232 232
233 list.push_back(particle1); 233 list.push_back(particle1);
234 list.push_back(particle2); 234 list.push_back(particle2);
235 list.push_back(kaon); 235 list.push_back(kaon);
236 236
237 PhaseSpaceGenerator::generateBiased(sqrtS, 237 PhaseSpaceGenerator::generateBiased(sqrtS, list, list.size()-3, angularSlope);
>> 238
>> 239 // INCL_DEBUG("NN Missing " << (kaon->getMomentum().theta()) * 180. / G4INCL::Math::pi << '\n');
238 240
239 fs->addModifiedParticle(particle1); 241 fs->addModifiedParticle(particle1);
240 fs->addModifiedParticle(particle2); 242 fs->addModifiedParticle(particle2);
241 fs->addCreatedParticle(kaon); 243 fs->addCreatedParticle(kaon);
242 for(Int_t i=0; i<nbr_pions; i++) fs->addCr 244 for(Int_t i=0; i<nbr_pions; i++) fs->addCreatedParticle(list[i]);
243 245
244 } 246 }
245 } 247 }
246 248