| 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 "G4INCLNpiToMissingStrangenessChannel 38 #include "G4INCLNpiToMissingStrangenessChannel.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 NpiToMissingStrangenessChanne 49 const G4double NpiToMissingStrangenessChannel::angularSlope = 1.;
50 50
51 NpiToMissingStrangenessChannel::NpiToMissing 51 NpiToMissingStrangenessChannel::NpiToMissingStrangenessChannel(Particle *p1, Particle *p2)
52 : particle1(p1), particle2(p2) 52 : particle1(p1), particle2(p2)
53 {} 53 {}
54 54
55 NpiToMissingStrangenessChannel::~NpiToMissin 55 NpiToMissingStrangenessChannel::~NpiToMissingStrangenessChannel(){}
56 56
57 void NpiToMissingStrangenessChannel::fillFin 57 void NpiToMissingStrangenessChannel::fillFinalState(FinalState *fs) {
58 58
59 const G4double sqrtS = KinematicsUtils::to 59 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(particle1, particle2); // MeV /!\/!\/!\.
60 // assert(sqrtS > 2.240); // ! > 2.109 Not sup 60 // assert(sqrtS > 2.240); // ! > 2.109 Not supposed to be under 2.244 GeV.
61 61
62 const G4int iso = ParticleTable::getIsospi 62 const G4int iso = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
63 // assert(iso == -3 || iso == -1 || iso == 1 | 63 // assert(iso == -3 || iso == -1 || iso == 1 || iso == 3);
64 G4int iso_system = 0; 64 G4int iso_system = 0;
65 G4int available_iso = 0; 65 G4int available_iso = 0;
66 G4int nbr_pions = 0; 66 G4int nbr_pions = 0;
67 G4int min_pions = 0; 67 G4int min_pions = 0;
68 G4int max_pions = 0; 68 G4int max_pions = 0;
69 69
70 Particle *pion_initial; 70 Particle *pion_initial;
71 Particle *nucleon_initial; 71 Particle *nucleon_initial;
72 72
73 if(particle1->isPion()){ 73 if(particle1->isPion()){
74 pion_initial = particle1; 74 pion_initial = particle1;
75 nucleon_initial = particle2; 75 nucleon_initial = particle2;
76 } 76 }
77 else{ 77 else{
78 pion_initial = particle2; 78 pion_initial = particle2;
79 nucleon_initial = particle1; 79 nucleon_initial = particle1;
80 } 80 }
81 const G4double pLab = 0.001*KinematicsUtil 81 const G4double pLab = 0.001*KinematicsUtils::momentumInLab(pion_initial, nucleon_initial); // GeV /!\/!\/!\.
82 82
83 G4double rdm = Random::shoot(); 83 G4double rdm = Random::shoot();
84 84
85 //G4int nbr_particle = 2; << 85 G4int nbr_particle = 2;
86 86
87 if(rdm < 0.35){ 87 if(rdm < 0.35){
88 // Lambda-K chosen 88 // Lambda-K chosen
89 nucleon_initial->setType(Lambda); 89 nucleon_initial->setType(Lambda);
90 available_iso = 1; 90 available_iso = 1;
91 min_pions = 3; 91 min_pions = 3;
92 max_pions = G4int((sqrtS-ParticleTable:: 92 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Lambda)-ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus));
93 } 93 }
94 else if((iso == 0 && rdm < 0.55) || rdm < 94 else if((iso == 0 && rdm < 0.55) || rdm < 0.5){
95 // N-K-Kb chosen 95 // N-K-Kb chosen
96 //nbr_particle++; << 96 nbr_particle++;
97 available_iso = 3; 97 available_iso = 3;
98 min_pions = 1; 98 min_pions = 1;
99 max_pions = G4int((sqrtS-ParticleTable:: 99 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(Proton)-2.*ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus));
100 } 100 }
101 else{ 101 else{
102 // Sigma-K chosen 102 // Sigma-K chosen
103 available_iso = 3; 103 available_iso = 3;
104 min_pions = 3; 104 min_pions = 3;
105 max_pions = G4int((sqrtS-ParticleTable:: 105 max_pions = G4int((sqrtS-ParticleTable::getINCLMass(SigmaMinus)-ParticleTable::getINCLMass(KZero)-10.)/ParticleTable::getINCLMass(PiPlus));
106 } 106 }
107 // Gaussian noise + mean 107 // Gaussian noise + mean value nbr pions fonction energy (choice)
108 G4double intermediaire = min_pions + Rando 108 G4double intermediaire = min_pions + Random::gauss(2) + std::sqrt(pLab-2.2);
109 nbr_pions = std::min(max_pions,std::max(mi 109 nbr_pions = std::min(max_pions,std::max(min_pions,G4int(intermediaire )));
110 110
111 available_iso += nbr_pions*2; 111 available_iso += nbr_pions*2;
112 #ifdef INCLXX_IN_GEANT4_MODE << 112 nbr_particle += nbr_pions;
113 // Erase the parent resonance information <<
114 particle1->setParentResonancePDGCode(0); <<
115 particle1->setParentResonanceID(0); <<
116 particle2->setParentResonancePDGCode(0); <<
117 particle2->setParentResonanceID(0); <<
118 #endif <<
119 //nbr_particle += nbr_pions; <<
120 113
121 ParticleList list; 114 ParticleList list;
122 ParticleType PionType = PiZero; 115 ParticleType PionType = PiZero;
123 const ThreeVector &rcol1 = pion_initial->g 116 const ThreeVector &rcol1 = pion_initial->getPosition();
124 const ThreeVector zero; 117 const ThreeVector zero;
125 118
126 // (pip piz pim) (sp sz s 119 // (pip piz pim) (sp sz sm) (L S Kb)
127 //pip_p 0.63 0.26 0.11 0.73 0.25 0. 120 //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. 121 //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); 122 G4bool pip_p = (std::abs(iso) == 3);
130 //piz_p 0.32 0.45 0.23 0.52 0.40 0. 123 //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( 124 G4bool piz_p = (ParticleTable::getIsospin(pion_initial->getType()) == 0);
132 //pim_p 0.18 0.37 0.45 0.20 0.63 0. 125 //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); 126 G4bool pim_p = (!pip_p && !piz_p);
134 127
135 for(Int_t i=0; i<nbr_pions; i++){ 128 for(Int_t i=0; i<nbr_pions; i++){
136 Particle *pion = new Particle(PionType,z 129 Particle *pion = new Particle(PionType,zero,rcol1);
137 if(available_iso-std::abs(iso-iso_system 130 if(available_iso-std::abs(iso-iso_system) >= 4){
138 rdm = Random::shoot(); 131 rdm = Random::shoot();
139 if((pip_p && rdm < 0.54) || (piz_p && 132 if((pip_p && rdm < 0.54) || (piz_p && rdm < 0.32) || (pim_p && rdm < 0.45)){
140 pion->setType(ParticleTable::getPion 133 pion->setType(ParticleTable::getPionType(G4int(Math::sign(iso))*2)); //pip/pip/pim
141 iso_system += 2*G4int(Math::sign(iso 134 iso_system += 2*G4int(Math::sign(iso));
142 available_iso -= 2; 135 available_iso -= 2;
143 } 136 }
144 else if((pip_p && rdm < 0.80) || (piz_ 137 else if((pip_p && rdm < 0.80) || (piz_p && rdm < 0.77) || (pim_p && rdm < 0.82)){
145 pion->setType(PiZero); 138 pion->setType(PiZero);
146 available_iso -= 2; 139 available_iso -= 2;
147 } 140 }
148 else{ 141 else{
149 pion->setType(ParticleTable::getPion 142 pion->setType(ParticleTable::getPionType(-G4int(Math::sign(iso))*2));
150 iso_system -= 2*G4int(Math::sign(iso 143 iso_system -= 2*G4int(Math::sign(iso));
151 available_iso -= 2; 144 available_iso -= 2;
152 } 145 }
153 } 146 }
154 else if(available_iso-std::abs(iso-iso_s 147 else if(available_iso-std::abs(iso-iso_system) == 2){
155 rdm = Random::shoot(); 148 rdm = Random::shoot();
156 if((pip_p && rdm < 0.26/0.37 && (Math: 149 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: 150 (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: 151 (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); 152 pion->setType(PiZero);
160 available_iso -= 2; 153 available_iso -= 2;
161 } 154 }
162 else{ 155 else{
163 pion->setType(ParticleTable::getPion 156 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
164 iso_system += Math::sign(iso-iso_sys 157 iso_system += Math::sign(iso-iso_system)*2;
165 available_iso -= 2; 158 available_iso -= 2;
166 } 159 }
167 } 160 }
168 else if(available_iso-std::abs(iso-iso_s 161 else if(available_iso-std::abs(iso-iso_system) == 0){
169 pion->setType(ParticleTable::getPionTy 162 pion->setType(ParticleTable::getPionType(Math::sign(iso-iso_system)*2));
170 iso_system += Math::sign(iso-iso_syste 163 iso_system += Math::sign(iso-iso_system)*2;
171 available_iso -= 2; 164 available_iso -= 2;
172 } 165 }
173 else INCL_ERROR("Pion Generation Problem 166 else INCL_ERROR("Pion Generation Problem in NpiToMissingStrangenessChannel" << '\n');
174 list.push_back(pion); 167 list.push_back(pion);
175 } 168 }
176 169
177 if(nucleon_initial->isLambda()){ // K-Lamb 170 if(nucleon_initial->isLambda()){ // K-Lambda
178 // assert(available_iso == 1); 171 // assert(available_iso == 1);
179 pion_initial->setType(ParticleTable::get 172 pion_initial->setType(ParticleTable::getKaonType(iso-iso_system));
180 } 173 }
181 else if(min_pions == 1){ // N-K-Kb chosen 174 else if(min_pions == 1){ // N-K-Kb chosen
182 // assert(available_iso == 3); 175 // assert(available_iso == 3);
183 Particle *antikaon = new Particle(KMinus 176 Particle *antikaon = new Particle(KMinus,zero,rcol1);
184 if(std::abs(iso-iso_system) == 3){ 177 if(std::abs(iso-iso_system) == 3){
185 pion_initial->setType(ParticleTable::g 178 pion_initial->setType(ParticleTable::getKaonType((iso-iso_system)/3));
186 nucleon_initial->setType(ParticleTable 179 nucleon_initial->setType(ParticleTable::getNucleonType((iso-iso_system)/3));
187 antikaon->setType(ParticleTable::getAn 180 antikaon->setType(ParticleTable::getAntiKaonType((iso-iso_system)/3));
188 } 181 }
189 else if(std::abs(iso-iso_system) == 1){ 182 else if(std::abs(iso-iso_system) == 1){ // equi-repartition
190 rdm = G4int(Random::shoot()*3.)-1; 183 rdm = G4int(Random::shoot()*3.)-1;
191 nucleon_initial->setType(ParticleTable 184 nucleon_initial->setType(ParticleTable::getNucleonType((G4int(rdm+0.5)*2-1)*(iso_system-iso)));
192 pion_initial->setType(ParticleTable::g 185 pion_initial->setType(ParticleTable::getKaonType((std::abs(rdm*2)-1)*(iso-iso_system)));
193 antikaon->setType(ParticleTable::getAn 186 antikaon->setType(ParticleTable::getAntiKaonType((G4int(rdm-0.5)*2+1)*(iso-iso_system)));
194 } 187 }
195 else INCL_ERROR("Isospin non-conservatio 188 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
196 list.push_back(antikaon); 189 list.push_back(antikaon);
197 nbr_pions += 1; // need for addCreatedPa 190 nbr_pions += 1; // need for addCreatedParticle loop
198 } 191 }
199 else{// Sigma-K 192 else{// Sigma-K
200 // assert(available_iso == 3); 193 // assert(available_iso == 3);
201 if(std::abs(iso-iso_system) == 3){ 194 if(std::abs(iso-iso_system) == 3){
202 pion_initial->setType(ParticleTable::g 195 pion_initial->setType(ParticleTable::getKaonType((iso-iso_system)/3));
203 nucleon_initial->setType(ParticleTable 196 nucleon_initial->setType(ParticleTable::getSigmaType((iso-iso_system)*2/3));
204 } 197 }
205 else if(std::abs(iso-iso_system) == 1){ 198 else if(std::abs(iso-iso_system) == 1){
206 rdm = Random::shoot(); 199 rdm = Random::shoot();
207 if((pip_p && rdm < 0.73) || (piz_p && 200 if((pip_p && rdm < 0.73) || (piz_p && rdm < 0.32) || (pim_p && rdm < 0.45)){
208 nucleon_initial->setType(SigmaZero); 201 nucleon_initial->setType(SigmaZero);
209 pion_initial->setType(ParticleTable: 202 pion_initial->setType(ParticleTable::getKaonType(iso-iso_system));
210 } 203 }
211 else{ 204 else{
212 nucleon_initial->setType(ParticleTab 205 nucleon_initial->setType(ParticleTable::getSigmaType((iso-iso_system)*2));
213 pion_initial->setType(ParticleTable: 206 pion_initial->setType(ParticleTable::getKaonType(iso_system-iso));
214 } 207 }
215 } 208 }
216 else INCL_ERROR("Isospin non-conservatio 209 else INCL_ERROR("Isospin non-conservation in NNToMissingStrangenessChannel" << '\n');
217 } 210 }
218 211
219 list.push_back(pion_initial); 212 list.push_back(pion_initial);
220 list.push_back(nucleon_initial); 213 list.push_back(nucleon_initial);
221 214
222 PhaseSpaceGenerator::generateBiased(sqrtS, 215 PhaseSpaceGenerator::generateBiased(sqrtS, list, list.size()-1, angularSlope);
223 216
224 fs->addModifiedParticle(pion_initial); 217 fs->addModifiedParticle(pion_initial);
225 fs->addModifiedParticle(nucleon_initial); 218 fs->addModifiedParticle(nucleon_initial);
226 for(Int_t i=0; i<nbr_pions; i++) fs->addCr 219 for(Int_t i=0; i<nbr_pions; i++) fs->addCreatedParticle(list[i]);
227 220
228 } 221 }
229 } 222 }
230 223