source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4NucleiModel.cc @ 1315

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 46.4 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// $Id: G4NucleiModel.cc,v 1.48 2010/05/26 18:29:28 dennis Exp $
26// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
27//
28// 20100112  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
29// 20100114  M. Kelsey -- Use G4ThreeVector for position
30// 20100309  M. Kelsey -- Use new generateWithRandomAngles for theta,phi stuff;
31//              eliminate some unnecessary std::pow(), move ctor here
32// 20100407  M. Kelsey -- Replace std::vector<>::resize(0) with ::clear().
33//              Move output vectors from generateParticleFate() and
34//              ::generateInteractionPartners() to be data members; return via
35//              const-ref instead of by value.
36// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
37// 20100418  M. Kelsey -- Reference output particle lists via const-ref
38// 20100421  M. Kelsey -- Replace hardwired p/n masses with G4PartDef's
39// 20100517  M. Kelsey -- Use G4CascadeINterpolator for cross-section
40//              calculations.  use G4Cascade*Channel for total xsec in pi-N
41//              N-N channels.  Move absorptionCrossSection() from SpecialFunc.
42
43//#define CHC_CHECK
44
45#include "G4NucleiModel.hh"
46#include "G4CascadeInterpolator.hh"
47#include "G4CascadeNNChannel.hh"
48#include "G4CascadeNPChannel.hh"
49#include "G4CascadePPChannel.hh"
50#include "G4CascadePiMinusNChannel.hh"
51#include "G4CascadePiMinusPChannel.hh"
52#include "G4CascadePiPlusNChannel.hh"
53#include "G4CascadePiPlusPChannel.hh"
54#include "G4CascadePiZeroNChannel.hh"
55#include "G4CascadePiZeroPChannel.hh"
56#include "G4CollisionOutput.hh"
57#include "G4ElementaryParticleCollider.hh"
58#include "G4HadTmpUtil.hh"
59#include "G4InuclNuclei.hh"
60#include "G4InuclParticleNames.hh"
61#include "G4InuclSpecialFunctions.hh"
62#include "G4LorentzConvertor.hh"
63#include "G4Neutron.hh"
64#include "G4NucleiProperties.hh"
65#include "G4Proton.hh"
66
67using namespace G4InuclParticleNames;
68using namespace G4InuclSpecialFunctions;
69
70
71typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
72
73G4NucleiModel::G4NucleiModel() : verboseLevel(0) {
74  if (verboseLevel > 3) {
75    G4cout << " >>> G4NucleiModel::G4NucleiModel" << G4endl;
76  }
77}
78
79G4NucleiModel::G4NucleiModel(G4InuclNuclei* nuclei) : verboseLevel(0) {
80  generateModel(nuclei->getA(), nuclei->getZ());
81}
82
83
84void 
85G4NucleiModel::generateModel(G4double a, G4double z) {
86  if (verboseLevel > 3) {
87    G4cout << " >>> G4NucleiModel::generateModel" << G4endl;
88  }
89
90  const G4double AU = 1.7234;
91  const G4double cuu = 3.3836;
92  //  const G4double convertToFermis = 2.8197;
93  const G4double pf_coeff = 1.932;
94  const G4double pion_vp = 0.007; // in GeV
95  const G4double pion_vp_small = 0.007; 
96  const G4double radForSmall = 8.0; // fermi
97  const G4double piTimes4thirds = 4.189; // 4 Pi/3
98  const G4double mproton = G4Proton::Definition()->GetPDGMass() / GeV;
99  const G4double mneutron = G4Neutron::Definition()->GetPDGMass() / GeV;
100  const G4double alfa3[3] = { 0.7, 0.3, 0.01 }; // listing zone radius
101  const G4double alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 };
102
103  A = a;
104  Z = z;
105  neutronNumber = a - z;
106  protonNumber = z;
107  neutronNumberCurrent = neutronNumber;
108  protonNumberCurrent = protonNumber;
109
110// Set binding energies
111//  G4double dm = bindingEnergy(a, z);
112  G4double dm = G4NucleiProperties::GetBindingEnergy(G4lrint(a), G4lrint(z));
113
114  binding_energies.push_back(0.001 * std::fabs(G4NucleiProperties::GetBindingEnergy(G4lrint(a-1), G4lrint(z-1)) - dm)); // for P
115  binding_energies.push_back(0.001 * std::fabs(G4NucleiProperties::GetBindingEnergy(G4lrint(a-1), G4lrint(z)) - dm)); // for N
116
117  G4double CU = cuu*G4cbrt(a); // half-density radius * 2.8197
118  G4double D1 = CU/AU;
119  G4double D = std::exp(-D1);   
120  G4double CU2 = 0.0; 
121
122  if (a > 4.5) {
123    std::vector<G4double> ur;
124    G4int icase = 0;
125
126    if (a > 99.5) {
127      number_of_zones = 6;
128      ur.push_back(-D1);
129
130      for (G4int i = 0; i < number_of_zones; i++) {
131        G4double y = std::log((1.0 + D) / alfa6[i] - 1.0);
132        zone_radii.push_back(CU + AU * y);
133        ur.push_back(y);
134      }
135
136    } else if (a > 11.5) {
137      number_of_zones = 3;
138      ur.push_back(-D1);
139
140      for (G4int i = 0; i < number_of_zones; i++) {
141        G4double y = std::log((1.0 + D)/alfa3[i] - 1.0);
142        zone_radii.push_back(CU + AU * y);
143        ur.push_back(y);
144      }
145
146    } else {
147      number_of_zones = 3;
148      icase = 1;
149      ur.push_back(0.0);
150 
151      G4double CU1 = CU * CU;
152      CU2 = std::sqrt(CU1 * (1.0 - 1.0 / a) + 6.4);
153
154      for (G4int i = 0; i < number_of_zones; i++) {
155        G4double y = std::sqrt(-std::log(alfa3[i]));
156        zone_radii.push_back(CU2 * y);
157        ur.push_back(y);
158      }
159    }
160
161    G4double tot_vol = 0.0;
162    std::vector<G4double> v;
163    std::vector<G4double> v1;
164
165    G4int i(0);
166    for (i = 0; i < number_of_zones; i++) {
167      G4double v0;
168
169      if (icase == 0) {
170        v0 = volNumInt(ur[i], ur[i + 1], CU, D1);
171
172      } else {
173        v0 = volNumInt1(ur[i], ur[i + 1], CU2);
174      };
175 
176      v.push_back(v0);
177      tot_vol += v0;
178
179      v0 = zone_radii[i]*zone_radii[i]*zone_radii[i];
180      if (i > 0) v0 -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
181
182      v1.push_back(v0);
183    }
184
185    // Protons
186    G4double dd0 = z/tot_vol/piTimes4thirds;
187    std::vector<G4double> rod;
188    std::vector<G4double> pf;
189    std::vector<G4double> vz;
190
191    for (i = 0; i < number_of_zones; i++) {
192      G4double rd = dd0 * v[i] / v1[i];
193      rod.push_back(rd);
194      G4double pff = pf_coeff * G4cbrt(rd);
195      pf.push_back(pff);
196      vz.push_back(0.5 * pff * pff / mproton + binding_energies[0]);
197    }
198
199    nucleon_densities.push_back(rod);
200    zone_potentials.push_back(vz);
201    fermi_momenta.push_back(pf);
202
203    // Neutrons
204    dd0 = (a - z)/tot_vol/piTimes4thirds;
205    rod.clear();
206    pf.clear();
207    vz.clear();
208
209    for (i = 0; i < number_of_zones; i++) {
210      G4double rd = dd0 * v[i] / v1[i];
211      rod.push_back(rd);
212      G4double pff = pf_coeff * G4cbrt(rd);
213      pf.push_back(pff);
214      vz.push_back(0.5 * pff * pff / mneutron + binding_energies[1]);
215    };
216
217    nucleon_densities.push_back(rod);
218    zone_potentials.push_back(vz);
219    fermi_momenta.push_back(pf);
220
221    // pion stuff (primitive)
222    std::vector<G4double> vp(number_of_zones, pion_vp);
223    zone_potentials.push_back(vp);
224
225    // kaon potential (primitive)
226    std::vector<G4double> kp(number_of_zones, -0.015);
227    zone_potentials.push_back(kp);
228
229    // hyperon potential (primitive)
230    std::vector<G4double> hp(number_of_zones, 0.03);
231    zone_potentials.push_back(hp);
232
233  } else { // a < 5
234    number_of_zones = 1;
235    G4double smallRad = radForSmall;
236    if (a == 4) smallRad *= 0.7;
237    zone_radii.push_back(smallRad);
238    G4double vol = 1.0 / piTimes4thirds / (zone_radii[0]*zone_radii[0]*zone_radii[0]);
239
240    // proton
241    std::vector<G4double> rod;
242    std::vector<G4double> pf;
243    std::vector<G4double> vz;
244    for (G4int i = 0; i < number_of_zones; i++) {
245      G4double rd = vol*z;
246      rod.push_back(rd);
247      G4double pff = pf_coeff * G4cbrt(rd);
248      pf.push_back(pff);
249      vz.push_back(0.5 * pff * pff / mproton + binding_energies[0]);
250    }
251
252    nucleon_densities.push_back(rod);
253    zone_potentials.push_back(vz);
254    fermi_momenta.push_back(pf);
255
256    // neutron
257    rod.clear();
258    pf.clear();
259    vz.clear();
260    for (G4int i = 0; i < number_of_zones; i++) {
261      G4double rd = vol*(a-z);
262      rod.push_back(rd);
263      G4double pff = pf_coeff * G4cbrt(rd);
264      pf.push_back(pff);
265      vz.push_back(0.5 * pff * pff / mneutron + binding_energies[1]);
266    }
267
268    nucleon_densities.push_back(rod);
269    zone_potentials.push_back(vz);
270    fermi_momenta.push_back(pf);
271
272    // pion (primitive)
273    std::vector<G4double> vp(number_of_zones, pion_vp_small);
274    zone_potentials.push_back(vp);
275 
276    // kaon potential (primitive)
277    std::vector<G4double> kp(number_of_zones, -0.015);
278    zone_potentials.push_back(kp);
279
280    // hyperon potential (primitive)
281    std::vector<G4double> hp(number_of_zones, 0.03);
282    zone_potentials.push_back(hp);
283  }
284
285  nuclei_radius = zone_radii[zone_radii.size() - 1];
286
287  /*
288  // Print nuclear radii and densities
289  G4cout << " For A = " << a << " zone radii = ";
290  for (G4int i = 0; i < number_of_zones; i++) G4cout << zone_radii[i] << "  ";
291  G4cout << "  " << G4endl;
292
293  G4cout << " proton densities: ";
294  for (G4int i = 0; i < number_of_zones; i++)
295     G4cout << nucleon_densities[0][i] << "  ";
296  G4cout << G4endl;
297
298  G4cout << " neutron densities: ";
299  for (G4int i = 0; i < number_of_zones; i++)
300     G4cout << nucleon_densities[1][i] << "  ";
301  G4cout << G4endl;
302
303  G4cout << " protons per shell " ;
304  G4double rinner = 0.0;
305  G4double router = 0.0;
306  G4double shellVolume = 0.0;
307  for (G4int i = 0; i < number_of_zones; i++) {  // loop over zones
308    router = zone_radii[i];
309    shellVolume = piTimes4thirds*(router*router*router - rinner*rinner*rinner);
310    G4cout << G4lrint(shellVolume*nucleon_densities[0][i]) << "  ";
311    rinner = router;
312  }
313  G4cout << G4endl;
314
315  G4cout << " neutrons per shell " ;
316  rinner = 0.0;
317  router = 0.0;
318  shellVolume = 0.0;
319  for (G4int i = 0; i < number_of_zones; i++) {  // loop over zones
320    router = zone_radii[i];
321    shellVolume = piTimes4thirds*(router*router*router - rinner*rinner*rinner);
322    G4cout << G4lrint(shellVolume*nucleon_densities[1][i]) << "  ";
323    rinner = router;
324  }
325  G4cout << G4endl;
326  */
327}
328
329
330G4double G4NucleiModel::getFermiKinetic(G4int ip, G4int izone) const {
331  G4double ekin = 0.0;
332 
333  if (ip < 3 && izone < number_of_zones) {      // ip for proton/neutron only
334    G4double pf = fermi_momenta[ip - 1][izone]; 
335    G4double mass = G4InuclElementaryParticle::getParticleMass(ip);
336   
337    ekin = std::sqrt(pf * pf + mass * mass) - mass;
338  } 
339  return ekin;
340}
341
342
343G4double
344G4NucleiModel::volNumInt(G4double r1, G4double r2, 
345                         G4double, G4double d1) const {
346
347  if (verboseLevel > 3) {
348    G4cout << " >>> G4NucleiModel::volNumInt" << G4endl;
349  }
350
351  const G4double au3 = 5.11864;
352  const G4double epsilon = 1.0e-3;
353  const G4int itry_max = 1000;
354  G4double d2 = 2.0 * d1;
355  G4double dr = r2 - r1;
356  G4double fi = 0.5 * (r1 * (r1 + d2) / (1.0 + std::exp(r1)) + r2 * (r2 + d2) / (1.0 + std::exp(r2)));
357  G4double fun1 = fi * dr;
358  G4double fun;
359  G4double jc = 1;
360  G4double dr1 = dr;
361  G4int itry = 0;
362
363  while (itry < itry_max) {
364    dr *= 0.5;
365    itry++;
366
367    G4double r = r1 - dr;
368    fi = 0.0;
369    G4int jc1 = G4int(std::pow(2.0, jc - 1) + 0.1);
370
371    for (G4int i = 0; i < jc1; i++) { 
372      r += dr1; 
373      fi += r * (r + d2) / (1.0 + std::exp(r));
374    };
375
376    fun = 0.5 * fun1 + fi * dr;
377
378    if (std::fabs((fun - fun1) / fun) > epsilon) {
379      jc++;
380      dr1 = dr;
381      fun1 = fun;
382    } else {
383      break;
384    }
385
386  }
387
388  if (verboseLevel > 2){
389    if(itry == itry_max) G4cout << " volNumInt-> n iter " << itry_max << G4endl;
390  }
391
392  return au3 * (fun + d1 * d1 * std::log((1.0 + std::exp(-r1)) / (1.0 + std::exp(-r2))));
393}
394
395
396G4double
397G4NucleiModel::volNumInt1(G4double r1, G4double r2, 
398                          G4double cu2) const {
399  if (verboseLevel > 3) {
400    G4cout << " >>> G4NucleiModel::volNumInt1" << G4endl;
401  }
402
403  const G4double epsilon = 1.0e-3;
404  const G4int itry_max = 1000;
405
406  G4double dr = r2 - r1;
407  G4double fi = 0.5 * (r1 * r1 * std::exp(-r1 * r1) + r2 * r2 * std::exp(-r2 * r2));
408  G4double fun1 = fi * dr;
409  G4double fun;
410  G4double jc = 1;
411  G4double dr1 = dr;
412  G4int itry = 0;
413
414  while (itry < itry_max) {
415    dr *= 0.5;
416    itry++;
417    G4double r = r1 - dr;
418    fi = 0.0;
419    G4int jc1 = int(std::pow(2.0, jc - 1) + 0.1);
420
421    for (G4int i = 0; i < jc1; i++) { 
422      r += dr1; 
423      fi += r * r * std::exp(-r * r);
424    }
425
426    fun = 0.5 * fun1 + fi * dr; 
427
428    if (std::fabs((fun - fun1) / fun) > epsilon) {
429      jc++;
430      dr1 = dr;
431      fun1 = fun;
432
433    } else {
434      break;
435    }
436
437  }
438
439  if (verboseLevel > 2){
440    if (itry == itry_max) G4cout << " volNumInt1-> n iter " << itry_max << G4endl;
441  }
442
443  return cu2*cu2*cu2 * fun;
444}
445
446
447void G4NucleiModel::printModel() const {
448
449  if (verboseLevel > 3) {
450    G4cout << " >>> G4NucleiModel::printModel" << G4endl;
451  }
452
453  G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
454         << " proton binding energy " << binding_energies[0] << 
455    " neutron binding energy " << binding_energies[1] << G4endl
456         << " Nculei radius " << nuclei_radius << " number of zones " <<
457    number_of_zones << G4endl;
458
459  for (G4int i = 0; i < number_of_zones; i++)
460
461    G4cout << " zone " << i+1 << " radius " << zone_radii[i] << G4endl
462           << " protons: density " << getDensity(1,i) << " PF " << 
463      getFermiMomentum(1,i) << " VP " << getPotential(1,i) << G4endl
464           << " neutrons: density " << getDensity(2,i) << " PF " << 
465      getFermiMomentum(2,i) << " VP " << getPotential(2,i) << G4endl
466           << " pions: VP " << getPotential(3,i) << G4endl;
467}
468
469
470G4LorentzVector
471G4NucleiModel::generateNucleonMomentum(G4int type, G4int zone) const {
472  G4double pmod = getFermiMomentum(type, zone) * G4cbrt(inuclRndm());
473  G4double mass = G4InuclElementaryParticle::getParticleMass(type);
474
475  return generateWithRandomAngles(pmod, mass);
476}
477
478
479G4InuclElementaryParticle
480G4NucleiModel::generateNucleon(G4int type, G4int zone) const {
481  if (verboseLevel > 3) {
482    G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
483  }
484
485  G4LorentzVector mom = generateNucleonMomentum(type, zone);
486  return G4InuclElementaryParticle(mom, type);
487}
488
489
490G4InuclElementaryParticle
491G4NucleiModel::generateQuasiDeutron(G4int type1, G4int type2,
492                                    G4int zone) const {
493
494  if (verboseLevel > 3) {
495    G4cout << " >>> G4NucleiModel::generateQuasiDeutron" << G4endl;
496  }
497
498  // Quasideuteron consists of an unbound but associated nucleon pair
499 
500  // FIXME:  Why generate two separate nucleon momenta (randomly!) and
501  //         add them, instead of just throwing a net momentum for the
502  //         dinulceon state?  And why do I have to capture the two
503  //         return values into local variables?
504  G4LorentzVector mom1 = generateNucleonMomentum(type1, zone);
505  G4LorentzVector mom2 = generateNucleonMomentum(type2, zone);
506  G4LorentzVector dmom = mom1+mom2;
507
508  G4int dtype = 0;
509       if (type1*type2 == pro*pro) dtype = 111;
510  else if (type1*type2 == pro*neu) dtype = 112;
511  else if (type1*type2 == neu*neu) dtype = 122;
512
513  return G4InuclElementaryParticle(dmom, dtype);
514}
515
516
517void
518G4NucleiModel::generateInteractionPartners(G4CascadParticle& cparticle) {
519  if (verboseLevel > 3) {
520    G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
521  }
522
523  const G4double pi4by3 = 4.1887903; // 4 Pi / 3
524  const G4double small = 1.0e-10;
525  const G4double huge_num = 50.0;
526  const G4double pn_spec = 1.0;
527
528  //const G4double pn_spec = 0.5;
529
530  //const G4double young_cut = std::sqrt(10.0) * 0.5;
531  //const G4double young_cut = std::sqrt(10.0) * 0.45;
532  const G4double young_cut = std::sqrt(10.0) * 0.25;
533  //const G4double young_cut = std::sqrt(10.0) * 0.2;
534  //const G4double young_cut = std::sqrt(10.0) * 0.1;
535  //const G4double young_cut = 0.0;
536
537  thePartners.clear();          // Reset buffer for next cycle
538
539  G4int ptype = cparticle.getParticle().type();
540  G4int zone = cparticle.getCurrentZone();
541  G4double pmass = cparticle.getParticle().getMass();
542  G4LorentzVector pmom = cparticle.getParticle().getMomentum();
543  G4double r_in;
544  G4double r_out;
545
546  if (zone == number_of_zones) { // particle is outside
547    r_in = nuclei_radius;
548    r_out = 0.0;
549
550  } else if (zone == 0) { // particle is outside core
551    r_in = 0.0;
552    r_out = zone_radii[0];
553
554  } else {
555    r_in = zone_radii[zone - 1];
556    r_out = zone_radii[zone];
557  }; 
558
559  G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
560
561  if (verboseLevel > 2){
562    G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path << G4endl;
563  }
564
565  if (path < -small) { // something wrong
566    return;
567
568  } else if (std::fabs(path) < small) { // just on the boundary
569    path = 0.0; 
570
571    G4InuclElementaryParticle particle; // Dummy -- no type or momentum
572    thePartners.push_back(partner(particle, path));
573
574  } else { // normal case 
575    G4LorentzConvertor dummy_convertor;
576    dummy_convertor.setBullet(pmom, pmass);
577 
578    for (G4int ip = 1; ip < 3; ip++) { 
579      G4InuclElementaryParticle particle = generateNucleon(ip, zone);
580      dummy_convertor.setTarget(particle.getMomentum(), particle.getMass());
581      G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
582
583      // Total cross section converted from mb to fm**2
584      G4double csec = totalCrossSection(ekin, ptype * ip);
585
586      if(verboseLevel > 2){
587        G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl;
588      }
589
590      G4double dens = nucleon_densities[ip - 1][zone];
591      G4double rat = getRatio(ip);
592      G4double pw = -path * dens * csec * rat;
593
594      if (pw < -huge_num) pw = -huge_num;
595      pw = 1.0 - std::exp(pw);
596
597      if (verboseLevel > 2){
598        G4cout << " pw " << pw << " rat " << rat << G4endl;
599      }
600
601      G4double spath = path;
602
603      if (inuclRndm() < pw) {
604        spath = -1.0 / dens / csec / rat * std::log(1.0 - pw * inuclRndm());
605        if (cparticle.young(young_cut, spath)) spath = path;
606
607        if (verboseLevel > 2){
608          G4cout << " ip " << ip << " spath " << spath << G4endl;
609        }
610
611      };
612      if (spath < path) thePartners.push_back(partner(particle, spath));
613    }; 
614
615    if (verboseLevel > 2){
616      G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
617    }
618
619    if (cparticle.getParticle().pion()) { // absorption possible
620      if (verboseLevel > 2) {
621        G4cout << " trying quasi-deuterons with bullet: ";
622        cparticle.getParticle().printParticle();
623      }
624
625      std::vector<G4InuclElementaryParticle> qdeutrons(3);
626      std::vector<G4double> acsecs(3);
627
628      G4double tot_abs_csec = 0.0;
629      G4double abs_sec;
630      G4double vol = zone_radii[zone]*zone_radii[zone]*zone_radii[zone];
631
632      if (zone > 0) vol -= zone_radii[zone-1]*zone_radii[zone-1]*zone_radii[zone-1];
633      vol *= pi4by3; 
634
635      G4double rat  = getRatio(1); 
636      G4double rat1 = getRatio(2); 
637
638      G4InuclElementaryParticle ppd = generateQuasiDeutron(1, 1, zone);
639
640      if (ptype == 7 || ptype == 5) {
641        dummy_convertor.setTarget(ppd.getMomentum(), ppd.getMass());
642
643        G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
644
645        if (verboseLevel > 2) {
646          G4cout << " ptype=" << ptype << " using pp target" << G4endl;
647          ppd.printParticle();
648        }
649
650        abs_sec = absorptionCrossSection(ekin, ptype);
651        abs_sec *= nucleon_densities[0][zone] * nucleon_densities[0][zone]*
652          rat * rat * vol; 
653
654      } else {
655        abs_sec = 0.0;
656      }; 
657
658      // abs_sec = 0.0;
659      tot_abs_csec += abs_sec;
660      acsecs.push_back(abs_sec);
661      qdeutrons.push_back(ppd);
662
663      G4InuclElementaryParticle npd = generateQuasiDeutron(1, 2, zone);
664
665      dummy_convertor.setTarget(npd.getMomentum(), npd.getMass());
666
667      G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
668
669      if (verboseLevel > 2) {
670        G4cout << " using np target" << G4endl;
671        npd.printParticle();
672      }
673
674      abs_sec = absorptionCrossSection(ekin, ptype); 
675      abs_sec *= pn_spec * nucleon_densities[0][zone] * nucleon_densities[1][zone] *
676        rat * rat1 * vol; 
677      tot_abs_csec += abs_sec;
678      acsecs.push_back(abs_sec);
679      qdeutrons.push_back(npd);
680
681      G4InuclElementaryParticle nnd = generateQuasiDeutron(2, 2, zone);
682
683      if (ptype == 7 || ptype == 3) {
684        dummy_convertor.setTarget(nnd.getMomentum(), nnd.getMass());
685
686        G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
687
688        if (verboseLevel > 2) {
689          G4cout << " ptype=" << ptype << " using nn target" << G4endl;
690          nnd.printParticle();
691        }
692
693        abs_sec = absorptionCrossSection(ekin, ptype); 
694        abs_sec *= nucleon_densities[1][zone] * nucleon_densities[1][zone] *
695          rat1 * rat1 * vol; 
696
697      } else {
698        abs_sec = 0.0;
699      }; 
700
701      // abs_sec = 0.0;
702      tot_abs_csec += abs_sec;
703      acsecs.push_back(abs_sec);
704      qdeutrons.push_back(nnd);
705
706      if (verboseLevel > 2){
707        G4cout << " rod1 " << acsecs[0] << " rod2 " << acsecs[1] 
708               << " rod3 " << acsecs[2] << G4endl;
709      }
710
711      if (tot_abs_csec > small) {
712     
713        G4double pw = -path * tot_abs_csec;
714
715        if (pw < -huge_num) pw = -huge_num;
716        pw = 1.0 - std::exp(pw);
717
718        if (verboseLevel > 2){
719          G4cout << " pw " << pw << G4endl;
720        }
721
722        G4double apath = path;
723
724        if (inuclRndm() < pw) 
725          apath = -1.0 / tot_abs_csec * std::log(1.0 - pw * inuclRndm());
726
727        if (cparticle.young(young_cut, apath)) apath = path; 
728
729        if(verboseLevel > 2){
730          G4cout << " apath " << apath << " path " << path << G4endl;
731        }
732
733        if (apath < path) { // chose the qdeutron
734
735          G4double sl = inuclRndm() * tot_abs_csec;
736          G4double as = 0.0;
737
738          for (G4int i = 0; i < 3; i++) {
739            as += acsecs[i];
740
741            if (sl < as) { 
742
743              if (verboseLevel > 2){
744                G4cout << " deut type " << i << G4endl; 
745              }
746
747              thePartners.push_back(partner(qdeutrons[i], apath));
748
749              break;
750            };
751          };
752        };   
753      };
754    }; 
755
756    if(verboseLevel > 2){
757      G4cout << " after deutrons " << thePartners.size() << G4endl;
758    }
759 
760    if (thePartners.size() > 1) {               // Sort list by path length
761      std::sort(thePartners.begin(), thePartners.end(), sortPartners);
762    }
763
764    G4InuclElementaryParticle particle;         // Dummy for end of list
765    thePartners.push_back(partner(particle, path));
766  }
767 
768  return;
769}
770
771
772const std::vector<G4CascadParticle>&
773G4NucleiModel::generateParticleFate(G4CascadParticle& cparticle,
774                                    G4ElementaryParticleCollider* theElementaryParticleCollider) {
775  if (verboseLevel > 3)
776    G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
777
778  outgoing_cparticles.clear();          // Clear return buffer for this event
779
780  generateInteractionPartners(cparticle);       // Fills "thePartners" data
781
782  if (thePartners.empty()) { // smth. is wrong -> needs special treatment
783    G4cout << " generateParticleFate-> can not be here " << G4endl;
784    return outgoing_cparticles;
785  }
786
787  G4int npart = thePartners.size();
788
789  if (npart == 1) { // cparticle is on the next zone entry
790    // need to go here if particle outside nucleus ?
791    //
792    cparticle.propagateAlongThePath(thePartners[0].second);
793    cparticle.incrementCurrentPath(thePartners[0].second);
794    boundaryTransition(cparticle);
795    outgoing_cparticles.push_back(cparticle);
796   
797    if (verboseLevel > 2){
798      G4cout << " next zone " << G4endl;
799      cparticle.print();
800    }
801   
802  } else { // there are possible interactions
803   
804    G4ThreeVector old_position = cparticle.getPosition();
805   
806    G4InuclElementaryParticle bullet = cparticle.getParticle();
807   
808    G4bool no_interaction = true;
809   
810    G4int zone = cparticle.getCurrentZone();
811   
812    G4CollisionOutput output;
813
814    for (G4int i = 0; i < npart - 1; i++) {
815      if (i > 0) cparticle.updatePosition(old_position); 
816     
817      G4InuclElementaryParticle target = thePartners[i].first; 
818     
819      if (verboseLevel > 2){
820        if (target.quasi_deutron()) 
821          G4cout << " try absorption: target " << target.type() << " bullet " <<
822            bullet.type() << G4endl;
823      }
824
825      output.reset();
826      theElementaryParticleCollider->collide(&bullet, &target, output);
827     
828      if (verboseLevel > 2) output.printCollisionOutput();
829
830      // Don't need to copy list, as "output" isn't changed again below
831      const std::vector<G4InuclElementaryParticle>& outgoing_particles = 
832        output.getOutgoingParticles();
833     
834      if (passFermi(outgoing_particles, zone)) { // interaction
835        cparticle.propagateAlongThePath(thePartners[i].second);
836        G4ThreeVector new_position = cparticle.getPosition();
837       
838        for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 
839          G4CascadParticle temp(outgoing_particles[ip], new_position, zone, 0.0, 0);
840          outgoing_cparticles.push_back(temp);
841        }
842       
843        no_interaction = false;
844        current_nucl1 = 0;
845        current_nucl2 = 0;
846#ifdef CHC_CHECK
847        G4double out_charge = 0.0;
848       
849        for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) 
850          out_charge += outgoing_particles[ip].getCharge();
851       
852        G4cout << " multiplicity " << outgoing_particles.size() <<
853          " bul type " << bullet.type() << " targ type " << target.type() << 
854          G4endl << " initial charge " << bullet.getCharge() + target.getCharge() 
855               << " out charge " << out_charge << G4endl; 
856#endif
857       
858        if (verboseLevel > 2){
859          G4cout << " partner type " << target.type() << G4endl;
860        }
861       
862        if (target.nucleon()) {
863          current_nucl1 = target.type();
864         
865        } else {
866          if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
867         
868          current_nucl1 = (target.type() - 100) / 10;
869          current_nucl2 = target.type() - 100 - 10 * current_nucl1;
870        }   
871       
872        if (current_nucl1 == 1) {
873          protonNumberCurrent -= 1.0;
874         
875        } else {
876          neutronNumberCurrent -= 1.0;
877        }; 
878       
879        if (current_nucl2 == 1) {
880          protonNumberCurrent -= 1.0;
881         
882        } else if(current_nucl2 == 2) {
883          neutronNumberCurrent -= 1.0;
884        };
885       
886        break;
887      }; 
888    }  // loop over partners
889   
890    if (no_interaction) { // still no interactions
891      cparticle.updatePosition(old_position); 
892      cparticle.propagateAlongThePath(thePartners[npart - 1].second);
893      cparticle.incrementCurrentPath(thePartners[npart - 1].second);
894      boundaryTransition(cparticle);
895      outgoing_cparticles.push_back(cparticle);
896    };
897  }; 
898
899  return outgoing_cparticles;
900}
901
902G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles, 
903                                G4int zone) {
904  if (verboseLevel > 3) {
905    G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
906  }
907
908  for (G4int i = 0; i < G4int(particles.size()); i++) {
909
910    if (particles[i].nucleon()) {
911
912      if (verboseLevel > 2){
913        G4cout << " type " << particles[i].type() << " p " << particles[i].getMomModule()
914               << " pf " << fermi_momenta[particles[i].type() - 1][zone] << G4endl;
915      }
916
917      if (particles[i].getMomModule() < fermi_momenta[particles[i].type() - 1][zone]) {
918
919        if (verboseLevel > 2) {
920          G4cout << " rejected by fermi: type " << particles[i].type() << 
921            " p " << particles[i].getMomModule() << G4endl;
922        }
923
924        return false;
925      };
926    };
927  };
928  return true; 
929}
930
931void G4NucleiModel::boundaryTransition(G4CascadParticle& cparticle) {
932
933  if (verboseLevel > 3) {
934    G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
935  }
936
937  G4int zone = cparticle.getCurrentZone();
938
939  if (cparticle.movingInsideNuclei() && zone == 0) {
940    G4cout << " boundaryTransition-> in zone 0 " << G4endl;
941
942  } else {
943    G4LorentzVector mom = cparticle.getMomentum();
944    G4ThreeVector pos = cparticle.getPosition();
945
946    G4int type = cparticle.getParticle().type();
947
948    G4double pr = pos.dot(mom.vect());
949    G4double r = pos.mag();
950
951    pr /= r;
952
953    G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
954
955    G4double dv = getPotential(type,zone) - getPotential(type, next_zone);
956    //    G4cout << "Potentials for type " << type << " = "
957    //           << getPotential(type,zone) << " , "
958    //     << getPotential(type,next_zone) << G4endl;
959
960    G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr;
961
962    G4double p1r;
963
964    if (verboseLevel > 2){
965      G4cout << " type " << type << " zone " << zone
966             << " next " << next_zone
967             << " qv " << qv << " dv " << dv << G4endl;
968    }
969
970    if(qv <= 0.0) { // reflection
971      p1r = -pr;
972      cparticle.incrementReflectionCounter();
973
974    } else { // transition
975      p1r = std::sqrt(qv);
976      if(pr < 0.0) p1r = -p1r;
977      cparticle.updateZone(next_zone);
978      cparticle.resetReflection();
979    };
980 
981    G4double prr = (p1r - pr) / r; 
982
983    mom.setVect(mom.vect() + pos*prr);
984    cparticle.updateParticleMomentum(mom);
985  }; 
986}
987
988G4bool G4NucleiModel::worthToPropagate(const G4CascadParticle& cparticle) const {
989
990  if (verboseLevel > 3) {
991    G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
992  }
993
994  const G4double cut_coeff = 2.0;
995
996  G4bool worth = true;
997
998  if (cparticle.reflectedNow()) {
999    G4int zone = cparticle.getCurrentZone();
1000
1001    G4int ip = cparticle.getParticle().type();
1002
1003    if (cparticle.getParticle().getKineticEnergy() < cut_coeff *   
1004       getFermiKinetic(ip, zone)) worth = false; 
1005
1006    if (verboseLevel > 3) {
1007      G4cout << "ekin=" << cparticle.getParticle().getKineticEnergy()
1008             << " fermiKin=" << getFermiKinetic(ip, zone) << " : worth? "
1009             << worth << G4endl;
1010    }
1011  };
1012
1013  return worth;
1014}
1015
1016G4double G4NucleiModel::getRatio(G4int ip) const {
1017
1018  if (verboseLevel > 3) {
1019    G4cout << " >>> G4NucleiModel::getRatio" << G4endl;
1020  }
1021
1022  G4double rat;
1023  //  G4double ratm;
1024
1025  // Calculate number of protons and neutrons in local region
1026  //  G4double Athird = G4cbrt(A);
1027  //  G4double Nneut = Athird*(A-Z)/A;
1028  //  G4double Nprot = Athird*Z/A;
1029
1030  // Reduce number of
1031  if (ip == 1) {
1032    if (verboseLevel > 2){
1033      G4cout << " current " << protonNumberCurrent << " inp " << protonNumber << G4endl;
1034    }
1035
1036    rat = protonNumberCurrent/protonNumber;
1037
1038    // Calculate ratio modified for local region
1039    //    G4double deltaP = protonNumber - protonNumberCurrent;
1040    //    G4cout << " deltaP = " << deltaP << G4endl;
1041    //    ratm = std::max(0.0, (Nprot - deltaP)/Nprot);
1042
1043  } else {
1044    if (verboseLevel > 2){
1045      G4cout << " current " << neutronNumberCurrent << " inp " << neutronNumber << G4endl;
1046    }
1047
1048    rat = neutronNumberCurrent/neutronNumber;
1049
1050    // Calculate ratio modified for local region
1051    //    G4double deltaN = neutronNumber - neutronNumberCurrent;
1052    //   G4cout << " deltaN = " << deltaN << G4endl;
1053    //    ratm = std::max(0.0, (Nneut - deltaN)/Nneut);
1054  }
1055
1056  //  G4cout << " get ratio: ratm =  " << ratm << G4endl;
1057  return rat;
1058  //  return ratm;
1059}
1060
1061G4CascadParticle G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle) {
1062
1063  if (verboseLevel > 3) {
1064    G4cout << " >>> G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle)" << G4endl;
1065  }
1066
1067  const G4double large = 1000.0;
1068
1069  G4double s1 = std::sqrt(inuclRndm()); 
1070  G4double phi = randomPHI();
1071  G4double rz = nuclei_radius * s1;
1072
1073  G4ThreeVector pos(rz*std::cos(phi), rz*std::sin(phi),
1074                    -nuclei_radius*std::sqrt(1.0 - s1*s1));
1075 
1076  G4CascadParticle cpart(*particle, pos, number_of_zones, large, 0);
1077
1078  if (verboseLevel > 2) cpart.print();
1079
1080  return cpart;
1081}
1082
1083void G4NucleiModel::initializeCascad(G4InuclNuclei* bullet, 
1084                                     G4InuclNuclei* target,
1085                                     modelLists& output) {
1086
1087  if (verboseLevel > 3) {
1088    G4cout << " >>> G4NucleiModel::initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target)" << G4endl;
1089  }
1090
1091  const G4double large = 1000.0;
1092  const G4double max_a_for_cascad = 5.0;
1093  const G4double ekin_cut = 2.0;
1094  const G4double small_ekin = 1.0e-6;
1095  const G4double r_large2for3 = 62.0;
1096  const G4double r0forAeq3 = 3.92;
1097  const G4double s3max = 6.5;
1098  const G4double r_large2for4 = 69.14;
1099  const G4double r0forAeq4 = 4.16;
1100  const G4double s4max = 7.0;
1101  const G4int itry_max = 100;
1102
1103  // Convenient references to input buffer contents
1104  std::vector<G4CascadParticle>& casparticles = output.first;
1105  std::vector<G4InuclElementaryParticle>& particles = output.second;
1106
1107  casparticles.clear();         // Reset input buffer to fill this event
1108  particles.clear();
1109
1110  // first decide whether it will be cascad or compound final nuclei
1111  G4double ab = bullet->getA();
1112  G4double zb = bullet->getZ();
1113  G4double at = target->getA();
1114  G4double zt = target->getZ();
1115
1116  G4double massb = bullet->getMass();   // For creating LorentzVectors below
1117
1118  if (ab < max_a_for_cascad) {
1119
1120    G4double benb = 0.001 * G4NucleiProperties::GetBindingEnergy(G4lrint(ab), G4lrint(zb)) / ab;
1121    G4double bent = 0.001 * G4NucleiProperties::GetBindingEnergy(G4lrint(at), G4lrint(zt)) / at;
1122    G4double ben = benb < bent ? bent : benb;
1123
1124    if (bullet->getKineticEnergy()/ab > ekin_cut*ben) {
1125      G4int itryg = 0;
1126
1127      while (casparticles.size() == 0 && itryg < itry_max) {     
1128        itryg++;
1129
1130        if(itryg > 0) particles.clear();
1131     
1132        //    nucleons coordinates and momenta in nuclei rest frame
1133        std::vector<G4ThreeVector> coordinates;
1134        std::vector<G4LorentzVector> momentums;
1135     
1136        if (ab < 3.0) { // deutron, simplest case
1137          G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm());
1138          G4double s = 2.0 * inuclRndm() - 1.0;
1139          G4double r1 = r * std::sqrt(1.0 - s * s);
1140          G4double phi = randomPHI();
1141
1142          G4ThreeVector coord1(r1*std::cos(phi), r1*std::sin(phi), r*s);   
1143          coordinates.push_back(coord1);
1144
1145          coord1 *= -1.;
1146          coordinates.push_back(coord1);
1147
1148          G4double p = 0.0;
1149          G4bool bad = true;
1150          G4int itry = 0;
1151
1152          while (bad && itry < itry_max) {
1153            itry++;
1154            p = 456.0 * inuclRndm();
1155
1156            if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
1157               p * r > 312.0) bad = false;
1158          };
1159
1160          if (itry == itry_max)
1161            if (verboseLevel > 2){ 
1162              G4cout << " deutron bullet generation-> itry = " << itry_max << G4endl;   
1163            }
1164
1165          p = 0.0005 * p;
1166
1167          if (verboseLevel > 2){ 
1168            G4cout << " p nuc " << p << G4endl;
1169          }
1170
1171          G4LorentzVector mom = generateWithRandomAngles(p, massb);
1172
1173          momentums.push_back(mom);
1174          mom.setVect(-mom.vect());
1175          momentums.push_back(-mom);
1176        } else {
1177          G4int ia = int(ab + 0.5);
1178
1179          G4ThreeVector coord1;
1180
1181          G4bool badco = true;
1182
1183          G4int itry = 0;
1184       
1185          if (ab < 4.0) { // a == 3
1186            while (badco && itry < itry_max) {
1187              if (itry > 0) coordinates.clear();
1188              itry++;   
1189              G4int i(0);   
1190
1191              for (i = 0; i < 2; i++) {
1192                G4int itry1 = 0;
1193                G4double s, u, rho; 
1194                G4double fmax = std::exp(-0.5) / std::sqrt(0.5);
1195
1196                while (itry1 < itry_max) {
1197                  itry1++;
1198                  s = -std::log(inuclRndm());
1199                  u = fmax * inuclRndm();
1200                  rho = std::sqrt(s) * std::exp(-s);
1201
1202                  if (std::sqrt(s) * std::exp(-s) > u && s < s3max) {
1203                    s = r0forAeq3 * std::sqrt(s);
1204                    coord1 = generateWithRandomAngles(s).vect();
1205                    coordinates.push_back(coord1);
1206
1207                    if (verboseLevel > 2){
1208                      G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1209                    }
1210                    break;
1211                  };
1212                };
1213
1214                if (itry1 == itry_max) { // bad case
1215                  coord1.set(10000.,10000.,10000.);
1216                  coordinates.push_back(coord1);
1217                  break;
1218                };
1219              };
1220
1221              coord1 = -coordinates[0] - coordinates[1]; 
1222              if (verboseLevel > 2) {
1223                G4cout << " 3  r " << coord1.mag() << G4endl;
1224              }
1225
1226              coordinates.push_back(coord1);       
1227           
1228              G4bool large_dist = false;
1229
1230              for (i = 0; i < 2; i++) {
1231                for (G4int j = i+1; j < 3; j++) {
1232                  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1233
1234                  if (verboseLevel > 2) {
1235                    G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1236                  }
1237
1238                  if (r2 > r_large2for3) {
1239                    large_dist = true;
1240
1241                    break; 
1242                  };     
1243                };
1244
1245                if (large_dist) break;
1246              }; 
1247
1248              if(!large_dist) badco = false;
1249
1250            };
1251
1252          } else { // a >= 4
1253            G4double b = 3./(ab - 2.0);
1254            G4double b1 = 1.0 - b / 2.0;
1255            G4double u = b1 + std::sqrt(b1 * b1 + b);
1256            b = 1.0 / b;
1257            G4double fmax = (1.0 + u * b) * u * std::exp(-u);
1258         
1259            while (badco && itry < itry_max) {
1260
1261              if (itry > 0) coordinates.clear();
1262              itry++;
1263              G4int i(0);
1264           
1265              for (i = 0; i < ia-1; i++) {
1266                G4int itry1 = 0;
1267                G4double s, u; 
1268
1269                while (itry1 < itry_max) {
1270                  itry1++;
1271                  s = -std::log(inuclRndm());
1272                  u = fmax * inuclRndm();
1273
1274                  if (std::sqrt(s) * std::exp(-s) * (1.0 + b * s) > u && s < s4max) {
1275                    s = r0forAeq4 * std::sqrt(s);
1276                    coord1 = generateWithRandomAngles(s).vect();
1277                    coordinates.push_back(coord1);
1278
1279                    if (verboseLevel > 2) {
1280                      G4cout << " i " << i << " r " << coord1.mag() << G4endl;
1281                    }
1282
1283                    break;
1284                  };
1285                };
1286
1287                if (itry1 == itry_max) { // bad case
1288                  coord1.set(10000.,10000.,10000.);
1289                  coordinates.push_back(coord1);
1290                  break;
1291                };
1292              };
1293
1294              coord1 *= 0.0;    // Cheap way to reset
1295              for(G4int j = 0; j < ia -1; j++) coord1 -= coordinates[j];
1296
1297              coordinates.push_back(coord1);   
1298
1299              if (verboseLevel > 2){
1300                G4cout << " last r " << coord1.mag() << G4endl;
1301              }
1302           
1303              G4bool large_dist = false;
1304
1305              for (i = 0; i < ia-1; i++) {
1306                for (G4int j = i+1; j < ia; j++) {
1307             
1308                  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
1309
1310                  if (verboseLevel > 2){
1311                    G4cout << " i " << i << " j " << j << " r2 " << r2 << G4endl;
1312                  }
1313
1314                  if (r2 > r_large2for4) {
1315                    large_dist = true;
1316
1317                    break; 
1318                  };     
1319                };
1320
1321                if (large_dist) break;
1322              }; 
1323
1324              if (!large_dist) badco = false;
1325            };
1326          }; 
1327
1328          if(badco) {
1329            G4cout << " can not generate the nucleons coordinates for a "
1330                   << ab << G4endl;     
1331
1332            casparticles.clear();       // Return empty buffer on error
1333            particles.clear();
1334            return;
1335
1336          } else { // momentums
1337            G4double p, u, x;
1338            G4LorentzVector mom;
1339            //G4bool badp = True;
1340            G4int i(0);
1341
1342            for (i = 0; i < ia - 1; i++) {
1343              G4int itry = 0;
1344
1345              while(itry < itry_max) {
1346                itry++;
1347                u = -std::log(0.879853 - 0.8798502 * inuclRndm());
1348                x = u * std::exp(-u);
1349
1350                if(x > inuclRndm()) {
1351                  p = std::sqrt(0.01953 * u);
1352                  mom = generateWithRandomAngles(p, massb);
1353                  momentums.push_back(mom);
1354
1355                  break;
1356                };
1357              };
1358
1359              if(itry == itry_max) {
1360                G4cout << " can not generate proper momentum for a "
1361                       << ab << G4endl;
1362
1363                casparticles.clear();   // Return empty buffer on error
1364                particles.clear();
1365                return;
1366              }; 
1367
1368            };
1369            // last momentum
1370
1371            mom *= 0.;          // Cheap way to reset
1372            for(G4int j=0; j< ia-1; j++) mom -= momentums[j]; 
1373
1374            momentums.push_back(mom);
1375          }; 
1376        }
1377 
1378        // Coordinates and momenta at rest are generated, now back to the lab
1379        G4double rb = 0.0;
1380        G4int i(0);
1381
1382        for(i = 0; i < G4int(coordinates.size()); i++) {     
1383          G4double rp = coordinates[i].mag2();
1384
1385          if(rp > rb) rb = rp;
1386        };
1387
1388        // nuclei i.p. as a whole
1389        G4double s1 = std::sqrt(inuclRndm()); 
1390        G4double phi = randomPHI();
1391        G4double rz = (nuclei_radius + rb) * s1;
1392        G4ThreeVector global_pos(rz*std::cos(phi), rz*std::sin(phi),
1393                                 -(nuclei_radius+rb)*std::sqrt(1.0-s1*s1));
1394
1395        for (i = 0; i < G4int(coordinates.size()); i++) {
1396          coordinates[i] += global_pos;
1397        }; 
1398
1399        // all nucleons at rest
1400        std::vector<G4InuclElementaryParticle> raw_particles;
1401        G4int ia = G4int(ab + 0.5);
1402        G4int iz = G4int(zb + 0.5);
1403
1404        for (G4int ipa = 0; ipa < ia; ipa++) {
1405          G4int knd = ipa < iz ? 1 : 2;
1406          raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
1407        }; 
1408     
1409        G4InuclElementaryParticle dummy(small_ekin, 1);
1410        G4LorentzConvertor toTheBulletRestFrame;
1411        toTheBulletRestFrame.setBullet(dummy.getMomentum(), dummy.getMass());
1412        toTheBulletRestFrame.setTarget(bullet->getMomentum(),bullet->getMass());
1413        toTheBulletRestFrame.toTheTargetRestFrame();
1414
1415        particleIterator ipart;
1416
1417        for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
1418          ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum())); 
1419        };
1420
1421        // fill cascad particles and outgoing particles
1422
1423        for(G4int ip = 0; ip < G4int(raw_particles.size()); ip++) {
1424          G4LorentzVector mom = raw_particles[ip].getMomentum();
1425          G4double pmod = mom.rho();
1426          G4double t0 = -mom.vect().dot(coordinates[ip]) / pmod;
1427          G4double det = t0 * t0 + nuclei_radius * nuclei_radius
1428                       - coordinates[ip].mag2();
1429          G4double tr = -1.0;
1430
1431          if(det > 0.0) {
1432            G4double t1 = t0 + std::sqrt(det);
1433            G4double t2 = t0 - std::sqrt(det);
1434
1435            if(std::fabs(t1) <= std::fabs(t2)) {         
1436              if(t1 > 0.0) {
1437                if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1438              };
1439
1440              if(tr < 0.0 && t2 > 0.0) {
1441
1442                if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1443              };
1444
1445            } else {
1446              if(t2 > 0.0) {
1447
1448                if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
1449              };
1450
1451              if(tr < 0.0 && t1 > 0.0) {
1452                if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
1453              };
1454            }; 
1455
1456          };
1457
1458          if(tr >= 0.0) { // cascad particle
1459            coordinates[ip] += mom*tr / pmod;
1460            casparticles.push_back(G4CascadParticle(raw_particles[ip], 
1461                                                    coordinates[ip], 
1462                                                    number_of_zones, large, 0));
1463
1464          } else {
1465            particles.push_back(raw_particles[ip]); 
1466          }; 
1467        };
1468      };   
1469
1470      if(casparticles.size() == 0) {
1471        particles.clear();
1472
1473        G4cout << " can not generate proper distribution for " << itry_max
1474               << " steps " << G4endl;
1475      };   
1476    };
1477  };
1478
1479  if(verboseLevel > 2){
1480    G4int ip(0);
1481
1482    G4cout << " cascad particles: " << casparticles.size() << G4endl;
1483    for(ip = 0; ip < G4int(casparticles.size()); ip++) casparticles[ip].print();
1484
1485    G4cout << " outgoing particles: " << particles.size() << G4endl;
1486    for(ip = 0; ip < G4int(particles.size()); ip++)
1487      particles[ip].printParticle();
1488  }
1489
1490  return;       // Buffer has been filled
1491}
1492
1493
1494G4double G4NucleiModel::absorptionCrossSection(G4double ke, G4int type) const {
1495  if (type != pionPlus && type != pionMinus && type != pionZero) {
1496    G4cerr << " absorptionCrossSection only valid for incident pions" << G4endl;
1497    return 0.;
1498  }
1499
1500  // was 0.2 since the beginning, then changed to 1.0
1501  // now 0.1 to convert from mb to fm**2
1502  const G4double corr_fac = 1.0;
1503  G4double csec = 0.0;
1504 
1505  if (ke < 0.3) {
1506    csec = 0.1106 / std::sqrt(ke) - 0.8 + 0.08 / ((ke-0.123)*(ke-0.123) + 0.0056);
1507
1508  } else if (ke < 1.0) {
1509    csec = 3.6735 * (1.0-ke)*(1.0-ke);     
1510  };
1511
1512  if (csec < 0.0) csec = 0.0;
1513
1514  if (verboseLevel > 2) {
1515    G4cout << " ekin " << ke << " abs. csec " << corr_fac * csec << G4endl;   
1516  }
1517
1518  csec *= corr_fac;
1519
1520  return csec;
1521}
1522
1523G4double G4NucleiModel::totalCrossSection(G4double ke, G4int rtype) const
1524{
1525  static const G4double keScale[] = {
1526    0.0,  0.01, 0.013, 0.018, 0.024, 0.032, 0.042, 0.056, 0.075, 0.1,
1527    0.13, 0.18, 0.24,  0.32,  0.42,  0.56,  0.75,  1.0,   1.3,   1.8,
1528    2.4,  3.2,  4.2,   5.6,   7.5,  10.0,  13.0,  18.0,  24.0,  32.0};
1529  static const G4int NBINS = sizeof(keScale)/sizeof(G4double);
1530
1531  static G4CascadeInterpolator<NBINS> interp(keScale);
1532
1533  // Pion and nucleon scattering cross-sections are available elsewhere
1534  switch (rtype) {
1535  case pro*pro: return G4CascadePPChannel::getCrossSection(ke); break;
1536  case pro*neu: return G4CascadeNPChannel::getCrossSection(ke); break;
1537  case pip*pro: return G4CascadePiPlusPChannel::getCrossSection(ke); break;
1538  case neu*neu: return G4CascadeNNChannel::getCrossSection(ke); break;
1539  case pim*pro: return G4CascadePiMinusPChannel::getCrossSection(ke); break;
1540  case pip*neu: return G4CascadePiPlusNChannel::getCrossSection(ke); break;
1541  case pi0*pro: return G4CascadePiZeroPChannel::getCrossSection(ke); break;
1542  case pim*neu: return G4CascadePiMinusNChannel::getCrossSection(ke); break;
1543  case pi0*neu: return G4CascadePiZeroNChannel::getCrossSection(ke); break;
1544    // Remaining channels are handled locally until arrays are moved
1545  case kpl*pro:                     
1546  case k0*neu:  return interp.interpolate(ke, kpPtot); break;
1547  case kmi*pro:                     
1548  case k0b*neu: return interp.interpolate(ke, kmPtot); break;
1549  case kpl*neu:                     
1550  case k0*pro:  return interp.interpolate(ke, kpNtot); break;
1551  case kmi*neu:                     
1552  case k0b*pro: return interp.interpolate(ke, kmNtot); break;
1553  case lam*pro:                     
1554  case lam*neu:                     
1555  case s0*pro:                       
1556  case s0*neu:  return interp.interpolate(ke, lPtot); break;
1557  case sp*pro:                       
1558  case sm*neu:  return interp.interpolate(ke, spPtot); break;
1559  case sm*pro:                       
1560  case sp*neu:  return interp.interpolate(ke, smPtot); break;
1561  case xi0*pro:                     
1562  case xim*neu: return interp.interpolate(ke, xi0Ptot); break;
1563  case xim*pro:                     
1564  case xi0*neu: return interp.interpolate(ke, ximPtot); break;
1565  default:
1566    G4cout << " unknown collison type = " << rtype << G4endl; 
1567  }
1568
1569  return 0.;            // Failure
1570}
1571
1572// Initialize cross-section interpolation tables
1573
1574const G4double G4NucleiModel::kpPtot[30] = {
1575   10.0,  10.34, 10.44, 10.61, 10.82, 11.09, 11.43, 11.71, 11.75, 11.8,
1576   11.98, 12.28, 12.56, 12.48, 12.67, 14.48, 15.92, 17.83, 17.93, 17.88,
1577   17.46, 17.3,  17.3,  17.4,  17.4,  17.4,  17.4,  17.5,  17.7,  17.8};
1578
1579const G4double G4NucleiModel::kpNtot[30] = {
1580    6.64,  6.99,  7.09,  7.27,  7.48,  7.75,  8.1,  8.49,  8.84, 9.31,
1581    9.8,  10.62, 11.64, 13.08, 14.88, 16.60, 17.5, 18.68, 18.68, 18.29,
1582   17.81, 17.6,  17.6,  17.6,  17.6,  17.6,  17.7, 17.8,  17.9,  18.0};
1583
1584const G4double G4NucleiModel::kmPtot[30] = {
1585 1997.0, 1681.41, 1586.74, 1428.95, 1239.59, 987.12, 671.54, 377.85, 247.30, 75.54,
1586    71.08, 54.74,   44.08,   44.38,   45.45,  45.07,  41.04,  35.75,  33.22, 30.08,
1587    27.61, 26.5,    25.2,    24.0,    23.4,   22.8,   22.0,   21.3,   21.0,  20.9};
1588
1589const G4double G4NucleiModel::kmNtot[30] = {
1590    6.15,  6.93,  7.16,  7.55,  8.02,  8.65,  9.43, 10.36, 11.34, 12.64,
1591   14.01, 16.45, 19.32, 23.0,  27.6,  30.92, 29.78, 28.28, 25.62, 23.1,
1592   22.31, 21.9,  21.73, 21.94, 21.23, 20.5,  20.4,  20.2,  20.1,  20.0};
1593
1594const G4double G4NucleiModel::lPtot[30] = {
1595  300.0, 249.07, 233.8, 208.33, 177.78, 137.04, 86.11, 41.41, 28.86, 12.35,
1596   13.82, 16.76, 20.68,  25.9,   30.37,  31.56, 32.83, 34.5,  34.91, 35.11,
1597   35.03, 36.06, 35.13,  35.01,  35.0,   35.0,  35.0,  35.0,  35.0,  35.0};
1598
1599const G4double G4NucleiModel::spPtot[30] = {
1600  150.0, 146.0, 144.8, 142.8, 140.4, 137.2, 133.2, 127.6, 120.0, 110.0,
1601   98.06, 84.16, 72.28, 56.58, 43.22, 40.44, 36.14, 30.48, 31.53, 31.92,
1602   29.25, 28.37, 29.81, 33.15, 33.95, 34.0,  34.0,  34.0,  34.0,  34.0};
1603
1604const G4double G4NucleiModel::smPtot[30] = {
1605  937.0, 788.14, 743.48, 669.05, 579.74, 460.65, 311.79, 183.33, 153.65, 114.6,
1606  105.18, 89.54,  70.58,  45.5,   32.17,  32.54,  32.95,  33.49,  33.55,  33.87,
1607   34.02, 34.29,  33.93,  33.88,  34.0,   34.0,   34.0,   34.0,   34.0,   34.0};
1608
1609const G4double G4NucleiModel::xi0Ptot[30] = {
1610  16.0,  14.72, 14.34, 13.7,  12.93, 11.9,  10.62, 9.29, 8.3,   7.0,
1611   7.96,  9.56, 11.48, 14.04, 19.22, 25.29, 29.4, 34.8, 34.32, 33.33,
1612  31.89, 29.55, 27.89, 21.43, 17.0,  16.0,  16.0, 16.0, 16.0,  16.0};
1613
1614const G4double G4NucleiModel::ximPtot[30] = {
1615  33.0,  32.5,  32.35, 32.1,  31.8,  31.4,  30.9, 30.2, 29.25, 28.0,
1616  26.5,  24.6,  22.8,  20.78, 18.22, 19.95, 21.7, 24.0, 24.74, 25.95,
1617  27.59, 27.54, 23.16, 17.43, 12.94, 12.0,  12.0, 12.0, 12.0,  12.0};
Note: See TracBrowser for help on using the repository browser.