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

Last change on this file since 1196 was 1196, checked in by garnier, 15 years ago

update CVS release candidate geant4.9.3.01
File size: 38.0 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  G4CascadeMomentum mom;
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  G4CascadeMomentum mom = generateNucleon(type1, zone).getMomentum(); 
377
378  G4CascadeMomentum mom1 = generateNucleon(type2, zone).getMomentum();
379
380  G4CascadeMomentum dmom;
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  const G4CascadeMomentum& 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, 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    G4CascadeMomentum mom = cparticle.getMomentum();
826    const 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, 0);
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<G4CascadeMomentum> 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          G4CascadeMomentum mom;
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            G4CascadeMomentum mom;
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          G4CascadeMomentum 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          const G4CascadeMomentum& 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, 0));
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.