| 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 "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(Particle *p1, Particle *p2)
50 : particle1(p1), particle2(p2)
51 {}
52
53 NNbarToLLbarChannel::~NNbarToLLbarChannel(){}
54
55 void NNbarToLLbarChannel::fillFinalState(FinalState *fs) {
56 // this channel include all states with lambdas, sigmas and xis and their antiparticles
57
58 //brief ppbar
59 // p pbar -> l lbar (BFMM 121)
60 // ppbar -> l lbar pi0 (BFMM 113)
61 // ppbar -> splus pim lbar || sminusbar pim l (BFMM 136)
62 // ppbar -> sminus pip lbar || splusbar l pip (BFMM 146)
63 // ppbar -> sp spbar (BFMM 139)
64 // ppbar -> sm smbar (BFMM 149)
65 // ppbar -> szero szerobar (BFMM 144)
66 // ppbar -> ximinus ximinusbar (BFMM 101)
67 // ppbar -> szero lbar || szerobar l (BFMM 143)
68 //
69 //
70 //brief npbar
71 // n pbar -> l lbar pi- (BFMM 487)
72 // n pbar -> l sbarplus || lbar sminus (BFMM 488)
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 487)
81 // p nbar -> l sbarminus || lbar splus (same as BFMM 488)
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*KinematicsUtils::momentumInLab(particle1, particle2); //GeV
97 // ppbar cross sections
98
99 const std::vector<G4double> BFMM121 = {2.379, -2.738, -1.260, -1.915, 0.430, 1.437};
100 //const G4double Eth_PPbar_LLbar = 1.437;
101 const std::vector<G4double> BFMM113 = {-0.105, 0.000, -5.099, 0.188, -0.050, 1.820};
102 //const G4double Eth_PPbar_LLbar_pi0 = 1.820;
103 const std::vector<G4double> BFMM139 = {0.142, -0.291, -1.702, -0.058, 0.001, 1.851};
104 //const G4double Eth_PPbar_SpSpbar = 1.851;
105 const std::vector<G4double> BFMM149 = {1.855, -2.238, -1.002, -1.279, 0.252, 1.896};
106 //const G4double Eth_PPbar_SmSmbar = 1.896;
107 const std::vector<G4double> BFMM136 = {1.749, -2.506, -1.222, -1.262, 0.274, 2.042};
108 //const G4double Eth_PPbar_SpLbar_pim = 2.042;
109 const std::vector<G4double> BFMM146 = {1.037, -1.437, -1.155, -0.709, 0.138, 2.065};
110 //const G4double Eth_PPbar_SmLbar_pip = 2.065;
111 const std::vector<G4double> BFMM143 = {0.652, -1.006, -1.805, -0.537, 0.121, 1.653};
112
113
114
115 //const G4double Eth_PPbar_Szero_Lbar = 1.653;
116 //fixed due to limited data
117 G4double BFMM144;
118 if(plab > 2.0) BFMM144 = 0.008; //sigmazero sigmazerobar
119 else BFMM144 = 0.0;
120 G4double BFMM101;
121 if(plab > 2.8) BFMM101 = 0.002; //ximinus ximinusbar
122 else BFMM101 = 0.0;
123
124 // npbar cross sections (fixed due to limited data)
125 G4double BFMM487;
126 if(plab > 2.1) BFMM487 = 0.048; //llbar piminus
127 else BFMM487 = 0.0;
128 G4double BFMM488;
129 if(plab > 2.0) BFMM488 = 0.139; //lsigmaminus +cc
130 else BFMM488 = 0.0;
131
132 const G4double sqrtS = KinematicsUtils::totalEnergyInCM(nucleon, antinucleon);
133 const G4double totalppbar = KinematicsUtils::compute_xs(BFMM113, plab)
134 +KinematicsUtils::compute_xs(BFMM139, plab) +KinematicsUtils::compute_xs(BFMM136, plab)
135 +KinematicsUtils::compute_xs(BFMM146, plab)+KinematicsUtils::compute_xs(BFMM143, plab)
136 +KinematicsUtils::compute_xs(BFMM121, plab)+KinematicsUtils::compute_xs(BFMM149, plab)
137 +BFMM144 +BFMM101;
138 const G4double totalpnbar = BFMM487 + BFMM488;
139 const G4double rdm = Random::shoot();
140
141 G4bool thirdparticle = false; //set true if we have pion
142 ParticleType PionType;
143 //setting types of new particles
144 if(nucleon->getType()==Proton){
145 if(antinucleon->getType()==antiProton){ //ppbar case
146 if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)){ //llbar
147 nucleon->setType(Lambda);
148 antinucleon->setType(antiLambda);
149 }
150 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144){ //sigmazero sigmazerobar
151 nucleon->setType(SigmaZero);
152 antinucleon->setType(antiSigmaZero);
153 }
154 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101){ //ximinus ximinusbar
155 nucleon->setType(XiMinus);
156 antinucleon->setType(antiXiMinus);
157 }
158 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
159 +KinematicsUtils::compute_xs(BFMM113, plab)){ //llbar pi0
160 nucleon->setType(Lambda);
161 antinucleon->setType(antiLambda);
162 thirdparticle = true;
163 PionType = PiZero;
164 }
165 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
166 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)){ //splus lbar pim || sminusbar l pim
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 < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
182 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)
183 +KinematicsUtils::compute_xs(BFMM146, plab)){ //sminus lbar pip || splussbar l pip
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 < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
199 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)
200 +KinematicsUtils::compute_xs(BFMM146, plab)+KinematicsUtils::compute_xs(BFMM143, plab)){ //szero lbar || szerobar l
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 < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
212 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)
213 +KinematicsUtils::compute_xs(BFMM146, plab)+KinematicsUtils::compute_xs(BFMM143, plab)
214 +KinematicsUtils::compute_xs(BFMM139, plab)){ //sp spbar
215 nucleon->setType(SigmaPlus);
216 antinucleon->setType(antiSigmaPlus);
217 }
218 else if(rdm*totalppbar < KinematicsUtils::compute_xs(std::move(BFMM121), plab)+BFMM144+BFMM101
219 +KinematicsUtils::compute_xs(std::move(BFMM113), plab)+KinematicsUtils::compute_xs(std::move(BFMM136), plab)
220 +KinematicsUtils::compute_xs(std::move(BFMM146), plab)+KinematicsUtils::compute_xs(std::move(BFMM143), plab)
221 +KinematicsUtils::compute_xs(std::move(BFMM139), plab)+KinematicsUtils::compute_xs(std::move(BFMM149), plab)){ //sm smbar
222 nucleon->setType(SigmaMinus);
223 antinucleon->setType(antiSigmaMinus);
224 }
225 else{
226 INCL_ERROR("out of total ppbar sum in LLbar channel");
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(antiSigmaMinus); //charge +1
235 }
236 else{
237 nucleon->setType(antiLambda);
238 antinucleon->setType(SigmaPlus); //charge +1
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){ //nnbar case same as ppbar
251 if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)){ //llbar
252 nucleon->setType(Lambda);
253 antinucleon->setType(antiLambda);
254 }
255 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144){ //sigmazero sigmazerobar
256 nucleon->setType(SigmaZero);
257 antinucleon->setType(antiSigmaZero);
258 }
259 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101){ //ximinus ximinusbar
260 nucleon->setType(XiMinus);
261 antinucleon->setType(antiXiMinus);
262 }
263 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
264 +KinematicsUtils::compute_xs(BFMM113, plab)){ //llbar pi0
265 nucleon->setType(Lambda);
266 antinucleon->setType(antiLambda);
267 thirdparticle = true;
268 PionType = PiZero;
269 }
270 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
271 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)){ //splus lbar pim || sminusbar l pim
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 < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
287 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)
288 +KinematicsUtils::compute_xs(BFMM146, plab)){ //sminus lbar pip || splussbar l pip
289 G4double rdm2 = Random::shoot();
290 if(rdm2 > 0.5){
291 nucleon->setType(SigmaMinus); //charge -1
292 antinucleon->setType(antiLambda);
293 thirdparticle = true;
294 PionType = PiPlus;
295 }
296 else{
297 nucleon->setType(antiSigmaPlus); //charge -1
298 antinucleon->setType(Lambda);
299 thirdparticle = true;
300 PionType = PiPlus;
301 }
302 }
303 else if(rdm*totalppbar < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
304 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)
305 +KinematicsUtils::compute_xs(BFMM146, plab)+KinematicsUtils::compute_xs(BFMM143, plab)){ //szero lbar || szerobar l
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 < KinematicsUtils::compute_xs(BFMM121, plab)+BFMM144+BFMM101
317 +KinematicsUtils::compute_xs(BFMM113, plab) + KinematicsUtils::compute_xs(BFMM136, plab)
318 +KinematicsUtils::compute_xs(BFMM146, plab)+KinematicsUtils::compute_xs(BFMM143, plab)
319 +KinematicsUtils::compute_xs(BFMM139, plab)){ //sp spbar
320 nucleon->setType(SigmaPlus);
321 antinucleon->setType(antiSigmaPlus);
322 }
323 else if(rdm*totalppbar < KinematicsUtils::compute_xs(std::move(BFMM121), plab)+BFMM144+BFMM101
324 +KinematicsUtils::compute_xs(std::move(BFMM113), plab)+KinematicsUtils::compute_xs(std::move(BFMM136), plab)
325 +KinematicsUtils::compute_xs(std::move(BFMM146), plab)+KinematicsUtils::compute_xs(std::move(BFMM143), plab)
326 +KinematicsUtils::compute_xs(std::move(BFMM139), plab)+KinematicsUtils::compute_xs(std::move(BFMM149), plab)){ //sm smbar
327 nucleon->setType(SigmaMinus);
328 antinucleon->setType(antiSigmaMinus);
329 }
330 else{
331 INCL_ERROR("out of total nnbar sum in LLbar channel");
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); //charge -1
340 }
341 else{
342 nucleon->setType(antiLambda);
343 antinucleon->setType(SigmaMinus); //charge -1
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 particles
356
357 if(thirdparticle){ //three particles
358 ParticleList list;
359 list.push_back(nucleon);
360 list.push_back(antinucleon);
361 const ThreeVector &rcol = nucleon->getPosition();
362 const ThreeVector zero;
363 Particle *pion = new Particle(PionType,zero,rcol);
364 list.push_back(pion);
365
366 PhaseSpaceGenerator::generate(sqrtS, list);
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*sqrtS);
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::normVector(py);
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