source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4BertiniNucleiModel.cc@ 819

Last change on this file since 819 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 38.1 KB
Line 
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#include "G4BertiniNucleiModel.hh"
27#include "G4LorentzConvertor.hh"
28#include "G4CollisionOutput.hh"
29
30typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
31
32G4BertiniNucleiModel::G4BertiniNucleiModel()
33 : verboseLevel(2) {
34
35 if (verboseLevel > 3) {
36 G4cout << " >>> G4BertiniNucleiModel::G4BertiniNucleiModel" << G4endl;
37 }
38}
39
40void G4BertiniNucleiModel::generateModel(G4double a,
41 G4double z) {
42 verboseLevel = 2;
43
44 if (verboseLevel > 3) {
45 G4cout << " >>> G4BertiniNucleiModel::generateModel" << G4endl;
46 }
47
48 const G4double AU = 1.7234;
49 const G4double cuu = 3.3836;
50 const G4double one_third = 1.0 / 3.0;
51 const G4double oneBypiTimes4 = 0.0795775; /// 1 / 4 Pi
52 const G4double pf_coeff = 1.932;
53 const G4double pion_vp = 0.007;
54 const G4double pion_vp_small = 0.007;
55 const G4double radForSmall = 8.0; /// fermi
56 const G4double piTimes4thirds = 4.189; /// 4 Pi/3
57 const G4double mproton = 0.93827;
58 const G4double mneutron = 0.93957;
59 const G4double alfa3[3] = { 0.7, 0.3, 0.01 };
60 /// const G4double alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 }; /// Case of six layers.
61
62 A = a;
63 Z = z;
64 neutronNumber = a - z;
65 protonNumber = z;
66 neutronNumberCurrent = neutronNumber;
67 protonNumberCurrent = protonNumber;
68
69/// set binding energies
70
71 G4double dm = bindingEnergy(a, z);
72 binding_energies.push_back(0.001 * std::fabs(bindingEnergy(a - 1, z - 1) - dm)); /// for P
73 binding_energies.push_back(0.001 * std::fabs(bindingEnergy(a - 1, z) - dm)); /// for N
74
75 G4double CU = cuu * std::pow(a, one_third);
76 G4double D1 = CU / AU;
77 G4double D = std::exp(-D1);
78 G4double CU2 = 0.0;
79
80 if (a > 3.5) { /// a > 3
81 std::vector<G4double> ur;
82 G4int icase = 0;
83
84 if (a > 11.5) { /// a > 11
85 /// number_of_zones = 6;
86 number_of_zones = 3;
87 ur.push_back(-D1);
88
89 for (G4int i = 0; i < number_of_zones; i++) {
90 /// G4double y = std::log((1.0 + D) / alfa6[i] - 1.0);
91 G4double y = std::log((1.0 + D)/alfa3[i] - 1.0);
92 zone_radii.push_back(CU + AU * y);
93 ur.push_back(y);
94 };
95
96 } else {
97 number_of_zones = 3;
98 icase = 1;
99 ur.push_back(0.0);
100 G4double CU1 = CU * CU;
101 CU2 = std::sqrt(CU1 * (1.0 - 1.0 / a) + 6.4);
102
103 for (G4int i = 0; i < number_of_zones; i++) {
104 G4double y = std::sqrt(-std::log(alfa3[i]));
105 zone_radii.push_back(CU2 * y);
106 ur.push_back(y);
107 };
108 };
109
110 G4double tot_vol = 0.0;
111 std::vector<G4double> v;
112 std::vector<G4double> v1;
113 G4int i(0);
114
115 for (i = 0; i < number_of_zones; i++) {
116 G4double v0;
117
118 if (icase == 0) {
119 v0 = volNumInt(ur[i], ur[i + 1], CU, D1);
120
121 } else {
122 v0 = volNumInt1(ur[i], ur[i + 1], CU2);
123 };
124
125 v.push_back(v0);
126 tot_vol += v0;
127 v0 = (i == 0 ? std::pow(zone_radii[i], G4double(3)) : std::pow(zone_radii[i], G4double(3)) -
128 std::pow(zone_radii[i - 1], G4double(3)));
129 v1.push_back(v0);
130 };
131
132 /// proton
133 G4double dd0 = 3.0 * z * oneBypiTimes4 / tot_vol;
134 std::vector<G4double> rod;
135 std::vector<G4double> pf;
136 std::vector<G4double> vz;
137
138 for (i = 0; i < number_of_zones; i++) {
139 G4double rd = dd0 * v[i] / v1[i];
140 rod.push_back(rd);
141 G4double pff = pf_coeff * std::pow(rd, one_third);
142 pf.push_back(pff);
143 vz.push_back(0.5 * pff * pff / mproton + binding_energies[0]);
144 };
145
146 nucleon_densities.push_back(rod);
147 zone_potentials.push_back(vz);
148 fermi_momenta.push_back(pf);
149
150 /// neutron
151 dd0 = 3.0 * (a - z) * oneBypiTimes4 / tot_vol;
152 rod.resize(0);
153 pf.resize(0);
154 vz.resize(0);
155
156 for (i = 0; i < number_of_zones; i++) {
157 G4double rd = dd0 * v[i] / v1[i];
158 rod.push_back(rd);
159 G4double pff = pf_coeff * std::pow(rd, one_third);
160 pf.push_back(pff);
161 vz.push_back(0.5 * pff * pff / mneutron + binding_energies[1]);
162 };
163
164 nucleon_densities.push_back(rod);
165 zone_potentials.push_back(vz);
166 fermi_momenta.push_back(pf);
167 /// pion
168 std::vector<G4double> vp(number_of_zones, pion_vp);
169 zone_potentials.push_back(vp);
170
171 } else { /// a < 4
172 number_of_zones = 1;
173 zone_radii.push_back(radForSmall);
174 G4double vol = 1.0 / piTimes4thirds / std::pow(zone_radii[0], G4double(3));
175 std::vector<G4double> rod;
176 std::vector<G4double> pf;
177 std::vector<G4double> vz;
178 G4int i(0);
179
180 for (i = 0; i < number_of_zones; i++) {
181 G4double rd = vol;
182 rod.push_back(rd);
183 G4double pff = pf_coeff * std::pow(rd, one_third);
184 pf.push_back(pff);
185 vz.push_back(0.5 * pff * pff / mproton + binding_energies[0]);
186 };
187
188 nucleon_densities.push_back(rod);
189 zone_potentials.push_back(vz);
190 fermi_momenta.push_back(pf);
191
192 /// neutron
193 rod.resize(0);
194 pf.resize(0);
195 vz.resize(0);
196
197 for (i = 0; i < number_of_zones; i++) {
198 G4double rd = vol;
199 rod.push_back(rd);
200 G4double pff = pf_coeff * std::pow(rd, one_third);
201 pf.push_back(pff);
202 vz.push_back(0.5 * pff * pff / mneutron + binding_energies[1]);
203 };
204
205 nucleon_densities.push_back(rod);
206 zone_potentials.push_back(vz);
207 fermi_momenta.push_back(pf);
208
209 /// pion
210 std::vector<G4double> vp(number_of_zones, pion_vp_small);
211 zone_potentials.push_back(vp);
212 };
213
214 nuclei_radius = zone_radii[zone_radii.size() - 1];
215}
216
217G4double G4BertiniNucleiModel::volNumInt(G4double r1,
218 G4double r2,
219 G4double /*cu*/,
220 G4double d1) const {
221
222 if (verboseLevel > 3) {
223 G4cout << " >>> G4BertiniNucleiModel::volNumInt" << G4endl;
224 }
225
226 const G4double au3 = 5.11864;
227 const G4double epsilon = 1.0e-3;
228 const G4int itry_max = 1000;
229 G4double d2 = 2.0 * d1;
230 G4double dr = r2 - r1;
231 G4double fi = 0.5 * (r1 * (r1 + d2) / (1.0 + std::exp(r1)) + r2 * (r2 + d2) / (1.0 + std::exp(r2)));
232 G4double fun1 = fi * dr;
233 G4double fun;
234 G4double jc = 1;
235 G4double dr1 = dr;
236 G4int itry = 0;
237
238 while (itry < itry_max) {
239 dr *= 0.5;
240 itry++;
241 G4double r = r1 - dr;
242 fi = 0.0;
243 G4int jc1 = G4int(std::pow(G4double(2.0), jc - 1) + 0.1);
244
245 for (G4int i = 0; i < jc1; i++) {
246 r += dr1;
247 fi += r * (r + d2) / (1.0 + std::exp(r));
248 };
249
250 fun = 0.5 * fun1 + fi * dr;
251
252 if (std::fabs((fun - fun1) / fun) > epsilon) {
253 jc++;
254 dr1 = dr;
255 fun1 = fun;
256
257 } else {
258
259 break;
260
261 };
262 };
263
264 if (verboseLevel > 2){
265
266 if (itry == itry_max) G4cout << " volNumInt-> n iter " << itry_max << G4endl;
267
268 }
269
270 return au3 * (fun + d1 * d1 * std::log((1.0 + std::exp(-r1)) / (1.0 + std::exp(-r2))));
271}
272
273G4double G4BertiniNucleiModel::volNumInt1(G4double r1,
274 G4double r2,
275 G4double cu2) const {
276
277 if (verboseLevel > 3) {
278 G4cout << " >>> G4BertiniNucleiModel::volNumInt1" << G4endl;
279 }
280
281 const G4double epsilon = 1.0e-3;
282 const G4int itry_max = 1000;
283 G4double dr = r2 - r1;
284 G4double fi = 0.5 * (r1 * r1 * std::exp(-r1 * r1) + r2 * r2 * std::exp(-r2 * r2));
285 G4double fun1 = fi * dr;
286 G4double fun;
287 G4double jc = 1;
288 G4double dr1 = dr;
289 G4int itry = 0;
290
291 while (itry < itry_max) {
292 dr *= 0.5;
293 itry++;
294 G4double r = r1 - dr;
295 fi = 0.0;
296 G4int jc1 = G4int(std::pow(2.0, jc - 1) + 0.1);
297
298 for (G4int i = 0; i < jc1; i++) {
299 r += dr1;
300 fi += r * r * std::exp(-r * r);
301 };
302
303 fun = 0.5 * fun1 + fi * dr;
304
305 if (std::fabs((fun - fun1) / fun) > epsilon) {
306 jc++;
307 dr1 = dr;
308 fun1 = fun;
309
310 } else {
311
312 break;
313 };
314 };
315
316 if (verboseLevel > 2){
317 if (itry == itry_max) G4cout << " volNumInt1-> n iter " << itry_max << G4endl;
318
319 }
320
321 return std::pow(cu2, G4double(3)) * fun;
322}
323
324void G4BertiniNucleiModel::printModel() const {
325
326 if (verboseLevel > 3) {
327 G4cout << " >>> G4BertiniNucleiModel::printModel" << G4endl;
328 }
329
330 G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
331 << " proton binding energy " << binding_energies[0] <<
332 " neutron binding energy " << binding_energies[1] << G4endl
333 << " Nculei radius " << nuclei_radius << " number of zones " <<
334 number_of_zones << G4endl;
335
336 for (G4int i = 0; i < number_of_zones; i++)
337
338 G4cout << " zone " << i+1 << " radius " << zone_radii[i] << G4endl
339 << " protons: density " << getDensity(1,i) << " PF " <<
340 getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
341 << " neutrons: density " << getDensity(2,i) << " PF " <<
342 getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
343 << " pions: VP " << getPotential(3,i) << G4endl;
344
345}
346
347G4InuclElementaryParticle G4BertiniNucleiModel::generateNucleon(G4int type,
348 G4int zone) const {
349
350 if (verboseLevel > 3) {
351 G4cout << " >>> G4BertiniNucleiModel::generateNucleon" << G4endl;
352 }
353
354 const G4double one_third = 1.0 / 3.0;
355 ///G4double pmod = getFermiMomentum(type, zone) * std::pow(inuclRndm(), one_third);
356 G4double pmod = fermi_momenta[type - 1][zone] * std::pow(inuclRndm(), one_third);
357 std::vector<G4double> mom(4);
358 std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
359 G4double FI = randomPHI();
360 G4double pt = pmod * COS_SIN.second;
361 mom[1] = pt * std::cos(FI);
362 mom[2] = pt * std::sin(FI);
363 mom[3] = pmod * COS_SIN.first;
364
365 return G4InuclElementaryParticle(mom, type);
366}
367
368G4InuclElementaryParticle G4BertiniNucleiModel::generateQuasiDeutron(G4int type1,
369 G4int type2,
370 G4int zone) const {
371
372 if (verboseLevel > 3) {
373 G4cout << " >>> G4BertiniNucleiModel::generateQuasiDeutron" << G4endl;
374 }
375
376 std::vector<G4double> mom = generateNucleon(type1, zone).getMomentum();
377
378 std::vector<G4double> mom1 = generateNucleon(type2, zone).getMomentum();
379
380 std::vector<G4double> dmom(4);
381
382 for(G4int i = 1; i < 4; i++) dmom[i] = mom[i] + mom1[i];
383
384 G4int dtype = 0;
385
386 if (type1 * type2 == 1) {
387 dtype = 111;
388
389 } else if(type1 * type2 == 2) {
390
391 dtype = 112;
392
393 } else if(type1 * type2 == 4) {
394
395 dtype = 122;
396 };
397
398 // return G4InuclElementaryParticle(dmom, dtype);
399 return G4InuclElementaryParticle(dmom, dtype, 3);
400}
401
402partners G4BertiniNucleiModel::generateInteractionPartners(G4CascadParticle& cparticle) const {
403
404 if (verboseLevel > 3) {
405 G4cout << " >>> G4BertiniNucleiModel::generateInteractionPartners" << G4endl;
406 }
407
408 const G4double pi4by3 = 4.1887903; /// 4 Pi / 3
409 const G4double small = 1.0e-10;
410 const G4double huge_num = 50.0;
411 const G4double pn_spec = 1.0;
412 ///const G4double pn_spec = 0.5;
413 ///const G4double young_cut = std::sqrt(10.0) * 0.1;
414 ///const G4double young_cut = std::sqrt(10.0) * 0.5;
415 ///const G4double young_cut = std::sqrt(10.0) * 0.45;
416 const G4double young_cut = std::sqrt(10.0) * 0.25;
417 ///const G4double young_cut = std::sqrt(10.0) * 0.2;
418 ///const G4double young_cut = 0.0;
419
420 partners thePartners;
421
422 G4int ptype = cparticle.getParticle().type();
423 G4int zone = cparticle.getCurrentZone();
424 G4double pmass = cparticle.getParticle().getMass();
425 std::vector<G4double> pmom = cparticle.getParticle().getMomentum();
426 G4double r_in;
427 G4double r_out;
428
429 if (zone == number_of_zones) { /// particle is outside
430 r_in = nuclei_radius;
431 r_out = 0.0;
432
433 } else if (zone == 0) { /// particle is outside core
434 r_in = 0.0;
435 r_out = zone_radii[0];
436
437 } else {
438 r_in = zone_radii[zone - 1];
439 r_out = zone_radii[zone];
440 };
441
442 G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
443
444 if (verboseLevel > 2){
445 G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path << G4endl;
446 }
447
448 if (path < -small) { /// something wrong
449 return thePartners;
450
451 } else if(std::fabs(path) < small) { /// just on the bounday
452 path = 0.0;
453
454 G4InuclElementaryParticle particle;
455
456 particle.setModel(3);
457
458 thePartners.push_back(partner(particle, path));
459
460 } else { /// normal case
461 std::vector<G4InuclElementaryParticle> particles;
462 G4LorentzConvertor dummy_convertor;
463
464 dummy_convertor.setBullet(pmom, pmass);
465
466 for (G4int ip = 1; ip < 3; ip++) {
467 G4InuclElementaryParticle particle = generateNucleon(ip, zone);
468
469 dummy_convertor.setTarget(particle.getMomentum(), particle.getMass());
470
471 G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
472
473 G4double csec = crossSection(ekin, ptype * ip);
474
475 if (verboseLevel > 2){
476 G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl;
477 }
478
479 G4double dens = nucleon_densities[ip - 1][zone];
480 G4double rat = getRatio(ip);
481 /// double rat = 1.0;
482 G4double pw = -path * dens * csec * rat;
483
484 if (pw < -huge_num) pw = -huge_num;
485 pw = 1.0 - std::exp(pw);
486
487 if(verboseLevel > 2){
488 G4cout << " pw " << pw << " rat " << rat << G4endl;
489 }
490
491 G4double spath = path;
492
493 if (inuclRndm() < pw) {
494 spath = -1.0 / dens / csec / rat * std::log(1.0 - pw * inuclRndm());
495 if (cparticle.young(young_cut, spath)) spath = path;
496
497 if (verboseLevel > 2){
498 G4cout << " ip " << ip << " spath " << spath << G4endl;
499 }
500
501 };
502 if (spath < path) thePartners.push_back(partner(particle, spath));
503 };
504
505 if (verboseLevel > 2){
506 G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
507 }
508
509 if (cparticle.getParticle().pion()) { /// absorption possible
510 std::vector<G4InuclElementaryParticle> qdeutrons;
511 std::vector<G4double> acsecs;
512 G4double tot_abs_csec = 0.0;
513 G4double abs_sec;
514 G4double vol = std::pow(zone_radii[zone], G4double(3) );
515
516 if(zone > 0) vol -= std::pow(zone_radii[zone - 1], 3);
517 vol *= pi4by3;
518 G4double rat = getRatio(1);
519 G4double rat1 = getRatio(2);
520
521 G4InuclElementaryParticle ppd = generateQuasiDeutron(1, 1, zone);
522
523 if (ptype == 7 || ptype == 5) {
524 dummy_convertor.setTarget(ppd.getMomentum(), ppd.getMass());
525
526 G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
527
528 abs_sec = absorptionCrosSection(ekin, ptype);
529 abs_sec *= nucleon_densities[0][zone] * nucleon_densities[0][zone]*
530 rat * rat * vol;
531
532 } else {
533 abs_sec = 0.0;
534 };
535
536 /// abs_sec = 0.0;
537 tot_abs_csec += abs_sec;
538 acsecs.push_back(abs_sec);
539 qdeutrons.push_back(ppd);
540
541 G4InuclElementaryParticle npd = generateQuasiDeutron(1, 2, zone);
542
543 dummy_convertor.setTarget(npd.getMomentum(), npd.getMass());
544
545 G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
546
547 abs_sec = absorptionCrosSection(ekin, ptype);
548 abs_sec *= pn_spec * nucleon_densities[0][zone] * nucleon_densities[1][zone] *
549 rat * rat1 * vol;
550 tot_abs_csec += abs_sec;
551 acsecs.push_back(abs_sec);
552 qdeutrons.push_back(npd);
553
554 G4InuclElementaryParticle nnd = generateQuasiDeutron(2, 2, zone);
555
556 if (ptype == 7 || ptype == 3) {
557 dummy_convertor.setTarget(nnd.getMomentum(), nnd.getMass());
558
559 G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
560
561 abs_sec = absorptionCrosSection(ekin, ptype);
562 abs_sec *= nucleon_densities[1][zone] * nucleon_densities[1][zone] *
563 rat1 * rat1 * vol;
564
565 } else {
566 abs_sec = 0.0;
567 };
568
569 /// abs_sec = 0.0;
570 tot_abs_csec += abs_sec;
571 acsecs.push_back(abs_sec);
572 qdeutrons.push_back(nnd);
573
574 if (verboseLevel > 2){
575 G4cout << " rod1 " << acsecs[0] << " rod2 " << acsecs[1]
576 << " rod3 " << acsecs[2] << G4endl;
577 }
578
579 if (tot_abs_csec > small) {
580
581 G4double pw = -path * tot_abs_csec;
582
583 if (pw < -huge_num) pw = -huge_num;
584
585 pw = 1.0 - std::exp(pw);
586
587 if (verboseLevel > 2){
588 G4cout << " pw " << pw << G4endl;
589 }
590
591 G4double apath = path;
592
593 if (inuclRndm() < pw) apath = -1.0 / tot_abs_csec * std::log(1.0 - pw * inuclRndm());
594
595 if (cparticle.young(young_cut, apath)) apath = path;
596
597 if (verboseLevel > 2){
598 G4cout << " apath " << apath << " path " << path << G4endl;
599 }
600
601 if (apath < path) { /// chose the qdeutron
602 G4double sl = inuclRndm() * tot_abs_csec;
603 G4double as = 0.0;
604
605 for (G4int i = 0; i < 3; i++) {
606 as += acsecs[i];
607
608 if (sl < as) {
609
610 if (verboseLevel > 2){
611 G4cout << " deut type " << i << G4endl;
612 }
613
614 thePartners.push_back(partner(qdeutrons[i], apath));
615
616 break;
617 };
618 };
619 };
620 };
621 };
622
623 if (verboseLevel > 2){
624 G4cout << " after deutrons " << thePartners.size() << G4endl;
625 }
626
627 if (thePartners.size() > 1) { /// sort partners
628
629 for (G4int i = 0; i < G4int(thePartners.size()) - 1; i++) {
630 for(G4int j = i + 1; j < G4int(thePartners.size()); j++) {
631
632 if (thePartners[i].second > thePartners[j].second) {
633 G4InuclElementaryParticle particle = thePartners[i].first;
634
635 particle.setModel(3);
636 G4double pathi = thePartners[i].second;
637 thePartners[i] = partner(thePartners[j].first, thePartners[j].second);
638 thePartners[j] = partner(particle, pathi);
639 };
640 };
641 };
642 };
643
644 G4InuclElementaryParticle particle;
645
646 particle.setModel(3);
647
648 thePartners.push_back(partner(particle, path));
649 };
650
651 return thePartners;
652}
653
654std::vector<G4CascadParticle> G4BertiniNucleiModel::generateParticleFate(G4CascadParticle& cparticle,
655 G4ElementaryParticleCollider* theElementaryParticleCollider) {
656
657 if (verboseLevel > 3) {
658 G4cout << " >>> G4BertiniNucleiModel::generateParticleFate" << G4endl;
659 }
660
661 std::vector<G4CascadParticle> outgouing_cparticles;
662
663 partners thePartners = generateInteractionPartners(cparticle);
664
665 if(thePartners.empty()) { /// smth. is wrong -> needs special treatment
666
667 G4cout << " generateParticleFate-> can not be here " << G4endl;
668
669 } else {
670 G4int npart = thePartners.size();
671
672 if (npart == 1) { /// cparticle is on the next zone entry
673 cparticle.propagateAlongThePath(thePartners[0].second);
674 cparticle.incrementCurrentPath(thePartners[0].second);
675 boundaryTransition(cparticle);
676 outgouing_cparticles.push_back(cparticle);
677
678 if (verboseLevel > 2){
679 G4cout << " next zone " << G4endl;
680 cparticle.print();
681 }
682
683 } else { /// there are possible interactions
684 std::vector<G4double> old_position = cparticle.getPosition();
685 G4InuclElementaryParticle bullet = cparticle.getParticle();
686 G4bool no_interaction = true;
687 G4int zone = cparticle.getCurrentZone();
688
689 for (G4int i = 0; i < npart - 1; i++) {
690
691 if(i > 0) cparticle.updatePosition(old_position);
692
693 G4InuclElementaryParticle target = thePartners[i].first;
694
695 if (verboseLevel > 2){
696
697 if (target.quasi_deutron())
698 G4cout << " try absorption: target " << target.type() << " bullet " << bullet.type() << G4endl;
699 }
700
701 G4CollisionOutput output = theElementaryParticleCollider->collide(&bullet, &target);
702
703 if (verboseLevel > 2){
704 output.printCollisionOutput();
705 }
706
707 std::vector<G4InuclElementaryParticle> outgoing_particles =
708
709 output.getOutgoingParticles();
710 if(passFermi(outgoing_particles, zone)) { /// interaction
711 cparticle.propagateAlongThePath(thePartners[i].second);
712
713 std::vector<G4double> new_position = cparticle.getPosition();
714
715 for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++)
716 outgouing_cparticles.push_back(G4CascadParticle(outgoing_particles[ip],
717 new_position, zone, 0.0));
718
719 no_interaction = false;
720 current_nucl1 = 0;
721 current_nucl2 = 0;
722
723#ifdef CHC_CHECK
724 G4double out_charge = 0.0;
725
726 for (G4int ip = 0; ip < outgoing_particles.size(); ip++)
727 out_charge += outgoing_particles[ip].getCharge();
728
729 G4cout << " multiplicity " << outgoing_particles.size() <<
730 " bul type " << bullet.type() << " targ type " << target.type() <<
731 G4endl << " initial charge " << bullet.getCharge() + target.getCharge()
732 << " out charge " << out_charge << G4endl;
733#endif
734
735 if (verboseLevel > 2){
736 G4cout << " partner type " << target.type() << G4endl;
737 }
738
739 if (target.nucleon()) {
740 current_nucl1 = target.type();
741
742 } else {
743
744 if (verboseLevel > 2){
745 G4cout << " good absorption " << G4endl;
746 }
747
748 current_nucl1 = (target.type() - 100) / 10;
749 current_nucl2 = target.type() - 100 - 10 * current_nucl1;
750 };
751
752 if (current_nucl1 == 1) {
753 protonNumberCurrent -= 1.0;
754
755 } else {
756 neutronNumberCurrent -= 1.0;
757 };
758
759 if (current_nucl2 == 1) {
760 protonNumberCurrent -= 1.0;
761
762 } else if(current_nucl2 == 2) {
763 neutronNumberCurrent -= 1.0;
764 };
765
766 break;
767 };
768 };
769
770 if (no_interaction) { /// still now interactions
771 cparticle.updatePosition(old_position);
772 cparticle.propagateAlongThePath(thePartners[npart - 1].second);
773 cparticle.incrementCurrentPath(thePartners[npart - 1].second);
774 boundaryTransition(cparticle);
775 outgouing_cparticles.push_back(cparticle);
776 };
777 };
778 };
779
780 return outgouing_cparticles;
781}
782
783G4bool G4BertiniNucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles,
784 G4int zone) {
785 if (verboseLevel > 3) {
786 G4cout << " >>> G4BertiniNucleiModel::passFermi" << G4endl;
787 }
788
789 for (G4int i = 0; i < G4int(particles.size()); i++) {
790
791 if (particles[i].nucleon()) {
792
793 if (verboseLevel > 2){
794 G4cout << " type " << particles[i].type() << " p " << particles[i].getMomModule()
795 << " pf " << fermi_momenta[particles[i].type() - 1][zone] << G4endl;
796 }
797
798 if (particles[i].getMomModule() < fermi_momenta[particles[i].type() - 1][zone]) {
799
800 if (verboseLevel > 2) {
801 G4cout << " rejected by fermi: type " << particles[i].type() <<
802 " p " << particles[i].getMomModule() << G4endl;
803 }
804
805 return false;
806 };
807 };
808 };
809 return true;
810}
811
812void G4BertiniNucleiModel::boundaryTransition(G4CascadParticle& cparticle) {
813
814 if (verboseLevel > 3) {
815 G4cout << " >>> G4BertiniNucleiModel::boundaryTransition" << G4endl;
816 }
817
818 G4int zone = cparticle.getCurrentZone();
819
820 if (cparticle.movingInsideNuclei() && zone == 0) {
821 G4cout << " boundaryTransition-> in zone 0 " << G4endl;
822
823 } else {
824
825 std::vector<G4double> mom = cparticle.getMomentum();
826 std::vector<G4double> pos = cparticle.getPosition();
827 G4int type = cparticle.getParticle().type();
828 G4double pr = 0.0;
829 G4double r = 0.0;
830 G4int i(0);
831
832 for (i = 0; i < 3; i++) {
833 pr += pos[i] * mom[i + 1];
834 r += pos[i] * pos[i];
835 };
836
837 r = std::sqrt(r);
838 pr /= r;
839
840 G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
841 G4double dv = getPotential(type,zone) - getPotential(type, next_zone);
842 G4double qv = dv * dv - 2.0 * dv * mom[0] + pr * pr;
843 G4double p1r;
844
845 if (verboseLevel > 2){
846 G4cout << " type " << type << " zone " << zone
847 << " next " << next_zone << " qv " << qv
848 << " dv " << dv << G4endl;
849 }
850
851 if (qv <= 0.0) { /// reflection
852 p1r = -pr;
853 cparticle.incrementReflectionCounter();
854
855 } else { /// transition
856 p1r = std::sqrt(qv);
857 if(pr < 0.0) p1r = -p1r;
858 cparticle.updateZone(next_zone);
859 cparticle.resetReflection();
860 };
861
862 G4double prr = (p1r - pr) / r;
863
864 for (i = 0; i < 3; i++) mom[i + 1] += pos[i] * prr;
865 cparticle.updateParticleMomentum(mom);
866 };
867}
868
869G4bool G4BertiniNucleiModel::worthToPropagate(const G4CascadParticle& cparticle) const {
870
871 if (verboseLevel > 3) {
872 G4cout << " >>> G4BertiniNucleiModel::worthToPropagate" << G4endl;
873 }
874
875 const G4double cut_coeff = 2.0;
876 G4bool worth = true;
877
878 if (cparticle.reflectedNow()) {
879 G4int zone = cparticle.getCurrentZone();
880 G4int ip = cparticle.getParticle().type();
881
882 if (cparticle.getParticle().getKineticEnergy() < cut_coeff *
883 getFermiKinetic(ip, zone)) worth = false;
884 };
885
886 return worth;
887}
888
889G4double G4BertiniNucleiModel::getRatio(G4int ip) const {
890
891 if (verboseLevel > 3) {
892 G4cout << " >>> G4BertiniNucleiModel::getRatio" << G4endl;
893 }
894
895 G4double rat;
896
897 if (ip == 1) {
898 if (verboseLevel > 2){
899 G4cout << " current " << protonNumberCurrent << " inp " << protonNumber << G4endl;
900 }
901
902 rat = protonNumberCurrent / protonNumber;
903
904 } else {
905
906 if (verboseLevel > 2){
907 G4cout << " current " << neutronNumberCurrent << " inp " << neutronNumber << G4endl;
908 }
909 rat = neutronNumberCurrent / neutronNumber;
910 };
911
912 return rat;
913}
914
915G4CascadParticle G4BertiniNucleiModel::initializeCascad(G4InuclElementaryParticle* particle) {
916
917 if (verboseLevel > 3) {
918 G4cout << " >>> G4BertiniNucleiModel::initializeCascad(G4InuclElementaryParticle* particle)" << G4endl;
919 }
920
921 const G4double large = 1000.0;
922 G4double s1 = std::sqrt(inuclRndm());
923 G4double phi = randomPHI();
924 G4double rz = nuclei_radius * s1;
925 std::vector<G4double> pos(3);
926 pos[0] = rz * std::cos(phi);
927 pos[1] = rz * std::sin(phi);
928 pos[2] = -nuclei_radius * std::sqrt(1.0 - s1 * s1);
929
930 G4CascadParticle cpart(*particle, pos, number_of_zones, large);
931
932 if (verboseLevel > 2){
933 cpart.print();
934 }
935
936 return cpart;
937}
938
939std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> >
940G4BertiniNucleiModel::initializeCascad(G4InuclNuclei* bullet,
941 G4InuclNuclei* target) {
942
943 if (verboseLevel > 3) {
944 G4cout << " >>> G4BertiniNucleiModel::initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target)" << G4endl;
945 }
946
947 const G4double large = 1000.0;
948 const G4double max_a_for_cascad = 5.0;
949 const G4double ekin_cut = 2.0;
950 const G4double small_ekin = 1.0e-6;
951 const G4double r_large2for3 = 62.0;
952 const G4double r0forAeq3 = 3.92;
953 const G4double s3max = 6.5;
954 const G4double r_large2for4 = 69.14;
955 const G4double r0forAeq4 = 4.16;
956 const G4double s4max = 7.0;
957 const G4int itry_max = 100;
958
959 std::vector<G4CascadParticle> casparticles;
960 std::vector<G4InuclElementaryParticle> particles;
961
962
963 G4double ab = bullet->getA();
964 G4double zb = bullet->getZ();
965 G4double at = target->getA();
966 G4double zt = target->getZ();
967
968 if (ab < max_a_for_cascad) { /// first decide whether it will be cascad or compound final nuclei
969 G4double benb = 0.001 * bindingEnergy(ab, zb) / ab;
970 G4double bent = 0.001 * bindingEnergy(at, zt) / at;
971 G4double ben = benb < bent ? bent : benb;
972
973 if (bullet->getKineticEnergy() / ab > ekin_cut * ben) {
974 G4int itryg = 0;
975
976 while (casparticles.size() == 0 && itryg < itry_max) {
977 itryg++;
978
979 if (itryg > 0) particles.resize(0);
980
981
982 std::vector<std::vector<G4double> > coordinates; /// nucleons coordinates in nuclei rest frame
983 std::vector<std::vector<G4double> > momentums;
984
985 if (ab < 3.0) { /// deutron, simplest case
986 G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm());
987 G4double s = 2.0 * inuclRndm() - 1.0;
988 G4double r1 = r * std::sqrt(1.0 - s * s);
989 std::vector<G4double> coord1(3);
990 G4double phi = randomPHI();
991 coord1[0] = r1 * std::cos(phi);
992 coord1[1] = r1 * std::sin(phi);
993 coord1[2] = r * s;
994
995 coordinates.push_back(coord1);
996 G4int i(0);
997
998 for (i = 0; i < 3; i++) coord1[i] *= -1;
999 coordinates.push_back(coord1);
1000 G4double p = 0.0;
1001 G4bool bad = true;
1002 G4int itry = 0;
1003
1004 while (bad && itry < itry_max) {
1005 itry++;
1006 p = 456.0 * inuclRndm();
1007
1008 if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
1009 p * r > 312.0) bad = false;
1010 };
1011
1012 if (itry == itry_max)
1013
1014 if (verboseLevel > 2){
1015 G4cout << " deutron bullet generation-> itry = "
1016 << itry_max << G4endl;
1017 }
1018
1019 p = 0.0005 * p;
1020
1021 if (verboseLevel > 2){
1022 G4cout << " p nuc " << p << G4endl;
1023 }
1024
1025 std::vector<G4double> mom(4);
1026 std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
1027 G4double FI = randomPHI();
1028 G4double P1 = p * COS_SIN.second;
1029 mom[1] = P1 * std::cos(FI);
1030 mom[2] = P1 * std::sin(FI);
1031 mom[3] = p * COS_SIN.first;
1032
1033 momentums.push_back(mom);
1034
1035 for (i = 1; i < 4; i++) mom[i] *= -1;
1036 momentums.push_back(mom);
1037
1038 } else {
1039 G4int ia = G4int(ab + 0.5);
1040 std::vector<G4double> coord1(3);
1041 G4bool badco = true;
1042 G4int itry = 0;
1043
1044 if (ab < 4.0) { /// a == 3
1045
1046 while (badco && itry < itry_max) {
1047 if (itry > 0) coordinates.resize(0);
1048 itry++;
1049 G4int i(0);
1050
1051 for (i = 0; i < 2; i++) {
1052 G4int itry1 = 0;
1053 G4double s;
1054 G4double u;
1055 G4double rho;
1056 G4double fmax = std::exp(-0.5) / std::sqrt(0.5);
1057
1058 while (itry1 < itry_max) {
1059 itry1++;
1060 s = -std::log(inuclRndm());
1061 u = fmax * inuclRndm();
1062 rho = std::sqrt(s) * std::exp(-s);
1063
1064 if (std::sqrt(s) * std::exp(-s) > u && s < s3max) {
1065 s = r0forAeq3 * std::sqrt(s);
1066 std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
1067 u = s * COS_SIN.second;
1068 G4double phi = randomPHI();
1069 coord1[0] = u * std::cos(phi);
1070 coord1[1] = u * std::sin(phi);
1071 coord1[2] = s * COS_SIN.first;
1072 coordinates.push_back(coord1);
1073
1074 if (verboseLevel > 2){
1075 G4cout << " i " << i << " r " << std::sqrt(coord1[0] * coord1[0] +
1076 coord1[1] * coord1[1] +
1077 coord1[2] * coord1[2]) << G4endl;
1078 }
1079
1080 break;
1081 };
1082 };
1083
1084 if (itry1 == itry_max) { /// bad case
1085 coord1[0] = coord1[1] = coord1[2] = 10000.;
1086 coordinates.push_back(coord1);
1087
1088 break;
1089 };
1090 };
1091
1092 for (i = 0; i < 3; i++) coord1[i] = - coordinates[0][i] -
1093 coordinates[1][i];
1094
1095 if (verboseLevel > 2) {
1096 G4cout << " 3 r " << std::sqrt(coord1[0] * coord1[0] +
1097 coord1[1] * coord1[1] +
1098 coord1[2] * coord1[2]) << G4endl;
1099 }
1100
1101 coordinates.push_back(coord1);
1102
1103 G4bool large_dist = false;
1104
1105 for (i = 0; i < 2; i++) {
1106 for (G4int j = i+1; j < 3; j++) {
1107 G4double r2 = std::pow(coordinates[i][0] - coordinates[j][0], G4double(2)) +
1108 std::pow(coordinates[i][1] - coordinates[j][1], G4double(2)) +
1109 std::pow(coordinates[i][2] - coordinates[j][2], G4double(2));
1110
1111 if (verboseLevel > 2) {
1112 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1113 }
1114
1115 if (r2 > r_large2for3) {
1116 large_dist = true;
1117
1118 break;
1119 };
1120 };
1121
1122 if (large_dist) break;
1123 };
1124
1125 if (!large_dist) badco = false;
1126 };
1127
1128 } else { /// a >= 4
1129 G4double b = 3.0/(ab - 2.0);
1130 G4double b1 = 1.0 - b / 2.0;
1131 G4double u = b1 + std::sqrt(b1 * b1 + b);
1132 b = 1.0 / b;
1133 G4double fmax = (1.0 + u * b) * u * std::exp(-u);
1134
1135 while (badco && itry < itry_max) {
1136
1137 if (itry > 0) coordinates.resize(0);
1138 itry++;
1139 G4int i(0);
1140
1141 for (i = 0; i < ia-1; i++) {
1142 G4int itry1 = 0;
1143 G4double s;
1144 G4double u;
1145
1146 while (itry1 < itry_max) {
1147 itry1++;
1148 s = -std::log(inuclRndm());
1149 u = fmax * inuclRndm();
1150
1151 if (std::sqrt(s) * std::exp(-s) * (1.0 + b * s) > u && s < s4max) {
1152 s = r0forAeq4 * std::sqrt(s);
1153 std::pair<double, double> COS_SIN = randomCOS_SIN();
1154 u = s * COS_SIN.second;
1155 G4double phi = randomPHI();
1156 coord1[0] = u * std::cos(phi);
1157 coord1[1] = u * std::sin(phi);
1158 coord1[2] = s * COS_SIN.first;
1159 coordinates.push_back(coord1);
1160
1161 if (verboseLevel > 2) {
1162 G4cout << " i " << i << " r " << std::sqrt(coord1[0] * coord1[0] +
1163 coord1[1] * coord1[1] +
1164 coord1[2] * coord1[2]) << G4endl;
1165 }
1166 break;
1167 };
1168 };
1169
1170 if (itry1 == itry_max) { /// bad case
1171 coord1[0] = coord1[1] = coord1[2] = 10000.0;
1172 coordinates.push_back(coord1);
1173
1174 break;
1175 };
1176 };
1177
1178 for (i = 0; i < 3; i++) {
1179 coord1[i] = 0.0;
1180
1181 for (G4int j = 0; j < ia -1; j++) coord1[i] -= coordinates[j][i];
1182 };
1183
1184 coordinates.push_back(coord1);
1185
1186 if (verboseLevel > 2){
1187 G4cout << " last r " << std::sqrt(coord1[0] * coord1[0] +
1188 coord1[1] * coord1[1] +
1189 coord1[2] * coord1[2]) << G4endl;
1190 }
1191
1192 G4bool large_dist = false;
1193
1194 for (i = 0; i < ia-1; i++) {
1195
1196 for (G4int j = i+1; j < ia; j++) {
1197 G4double r2 = std::pow(coordinates[i][0] - coordinates[j][0], G4double(2)) +
1198 std::pow(coordinates[i][1]-coordinates[j][1], G4double(2)) +
1199 std::pow(coordinates[i][2] - coordinates[j][2], G4double(2));
1200
1201 if (verboseLevel > 2){
1202 G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1203 }
1204
1205 if (r2 > r_large2for4) {
1206 large_dist = true;
1207
1208 break;
1209 };
1210 };
1211
1212 if (large_dist) break;
1213 };
1214
1215 if (!large_dist) badco = false;
1216 };
1217 };
1218
1219 if (badco) {
1220 G4cout << " can not generate the nucleons coordinates for a " << ab <<
1221 G4endl;
1222
1223 return std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> >
1224 (casparticles, particles);
1225
1226 } else { /// momentums
1227 G4double p;
1228 G4double u;
1229 G4double x;
1230 std::vector<G4double> mom(4);
1231 /// G4bool badp = True;
1232 G4int i(0);
1233
1234 for (i = 0; i < ia - 1; i++) {
1235 G4int itry = 0;
1236
1237 while (itry < itry_max) {
1238 itry++;
1239 u = -std::log(0.879853 - 0.8798502 * inuclRndm());
1240 x = u * std::exp(-u);
1241
1242 if (x > inuclRndm()) {
1243 p = std::sqrt(0.01953 * u);
1244 std::pair<G4double, G4double> COS_SIN = randomCOS_SIN();
1245 G4double pt = p * COS_SIN.second;
1246 G4double phi = randomPHI();
1247 mom[1] = pt * std::cos(phi);
1248 mom[2] = pt * std::sin(phi);
1249 mom[3] = p * COS_SIN.first;
1250 momentums.push_back(mom);
1251
1252 break;
1253 };
1254 };
1255
1256 if (itry == itry_max) {
1257 G4cout << " can not generate proper momentum for a " << ab << G4endl;
1258
1259 return std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> >
1260 (casparticles, particles);
1261 };
1262 };
1263
1264 /// last momentum
1265 for (i = 1; i < 4; i++) {
1266 mom[i] = 0.;
1267
1268 for (G4int j = 0; j < ia -1; j++) mom[i] -= momentums[j][i];
1269 };
1270 momentums.push_back(mom);
1271 };
1272 };
1273
1274 /// coordinates and momentums at rest are generated, now back to the lab;
1275 G4double rb = 0.0;
1276 G4int i(0);
1277
1278 for (i = 0; i < G4int(coordinates.size()); i++) {
1279 G4double rp = std::sqrt(coordinates[i][0] * coordinates[i][0] +
1280 coordinates[i][1] * coordinates[i][1] +
1281 coordinates[i][2] * coordinates[i][2]);
1282 if (rp > rb) rb = rp;
1283 };
1284 /// nuclei i.p. as a whole
1285 G4double s1 = std::sqrt(inuclRndm());
1286 G4double phi = randomPHI();
1287 G4double rz = (nuclei_radius + rb) * s1;
1288 std::vector<double> global_pos(3);
1289 global_pos[0] = rz * std::cos(phi);
1290 global_pos[1] = rz * std::sin(phi);
1291 global_pos[2] = -(nuclei_radius + rb) * std::sqrt(1.0 - s1 * s1);
1292
1293 for (i = 0; i < G4int(coordinates.size()); i++) {
1294 coordinates[i][0] += global_pos[0];
1295 coordinates[i][1] += global_pos[1];
1296 coordinates[i][2] += global_pos[2];
1297 };
1298
1299 /// all nucleons at rest
1300 std::vector<G4InuclElementaryParticle> raw_particles;
1301 G4int ia = G4int(ab + 0.5);
1302 G4int iz = G4int(zb + 0.5);
1303
1304 for (G4int ipa = 0; ipa < ia; ipa++) {
1305 G4int knd = ipa < iz ? 1 : 2;
1306
1307 raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
1308 };
1309
1310 G4InuclElementaryParticle dummy(small_ekin, 1);
1311
1312 G4LorentzConvertor toTheBulletRestFrame;
1313
1314 toTheBulletRestFrame.setBullet(dummy.getMomentum(), dummy.getMass());
1315 toTheBulletRestFrame.setTarget(bullet->getMomentum(),bullet->getMass());
1316 toTheBulletRestFrame.toTheTargetRestFrame();
1317
1318 particleIterator ipart;
1319
1320 for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
1321 std::vector<G4double> mom =
1322 toTheBulletRestFrame.backToTheLab(ipart->getMomentum());
1323
1324 ipart->setMomentum(mom);
1325 };
1326
1327 /// fill cascad particles and outgoing particles
1328 for (G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
1329 std::vector<G4double> mom = raw_particles[ip].getMomentum();
1330 G4double pmod = std::sqrt(mom[1] * mom[1] + mom[2] * mom[2] + mom[3] * mom[3]);
1331 G4double t0 = -(mom[1] * coordinates[ip][0] + mom[2] * coordinates[ip][1] +
1332 mom[3] * coordinates[ip][2]) / pmod;
1333 G4double det = t0 * t0 + nuclei_radius * nuclei_radius -
1334 coordinates[ip][0] * coordinates[ip][0] -
1335 coordinates[ip][1] * coordinates[ip][1] -
1336 coordinates[ip][2] * coordinates[ip][2];
1337 G4double tr = -1.0;
1338
1339 if (det > 0.0) {
1340 G4double t1 = t0 + std::sqrt(det);
1341 G4double t2 = t0 - std::sqrt(det);
1342
1343 if (std::fabs(t1) <= std::fabs(t2)) {
1344
1345 if (t1 > 0.0) {
1346
1347 if (coordinates[ip][2] + mom[3] * t1 / pmod <= 0.0) tr = t1;
1348 };
1349
1350 if (tr < 0.0 && t2 > 0.0) {
1351
1352 if (coordinates[ip][2] + mom[3] * t2 / pmod <= 0.0) tr = t2;
1353 };
1354
1355 } else {
1356
1357 if (t2 > 0.0) {
1358
1359 if (coordinates[ip][2] + mom[3] * t2 / pmod <= 0.0) tr = t2;
1360
1361 };
1362
1363 if (tr < 0.0 && t1 > 0.0) {
1364
1365 if (coordinates[ip][2] + mom[3] * t1 / pmod <= 0.0) tr = t1;
1366 };
1367
1368 };
1369 };
1370
1371 if (tr >= 0.0) { /// cascad particle
1372 coordinates[ip][0] += mom[1] * tr / pmod;
1373 coordinates[ip][1] += mom[2] * tr / pmod;
1374 coordinates[ip][2] += mom[3] * tr / pmod;
1375 casparticles.push_back(G4CascadParticle(raw_particles[ip], coordinates[ip],
1376 number_of_zones, large));
1377 } else {
1378 particles.push_back(raw_particles[ip]);
1379 };
1380 };
1381 };
1382
1383 if (casparticles.size() == 0) {
1384 particles.resize(0);
1385
1386 G4cout << " can not generate proper distribution for " << itry_max << " steps " << G4endl;
1387 };
1388 };
1389 };
1390
1391 if (verboseLevel > 2){
1392 G4cout << " cascad particles: " << casparticles.size() << G4endl;
1393 G4int ip(0);
1394
1395 for (ip = 0; ip < G4int(casparticles.size()); ip++)
1396 casparticles[ip].print();
1397 G4cout << " outgoing particles: " << particles.size() << G4endl;
1398
1399 for (ip = 0; ip < G4int(particles.size()); ip++)
1400 particles[ip].printParticle();
1401 }
1402
1403 return std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> >
1404 (casparticles, particles);
1405}
Note: See TracBrowser for help on using the repository browser.