| 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