| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // J. M. Quesada (August 2008).
28 // Based on previous work by V. Lara
29 //
30 // Modified:
31 // 06.09.2008 JMQ Also external choice has been added for:
32 // - superimposed Coulomb barrier (if useSICB=true)
33 // 20.08.2010 V.Ivanchenko cleanup
34 //
35
36 #include "G4PreCompoundFragment.hh"
37 #include "G4KalbachCrossSection.hh"
38 #include "G4ChatterjeeCrossSection.hh"
39 #include "G4DeexPrecoParameters.hh"
40 #include "G4InterfaceToXS.hh"
41 #include "G4IsotopeList.hh"
42 #include "Randomize.hh"
43
44 G4PreCompoundFragment::G4PreCompoundFragment(const G4ParticleDefinition* p,
45 G4VCoulombBarrier* aCoulBarrier)
46 : G4VPreCompoundFragment(p, aCoulBarrier)
47 {}
48
49 G4double G4PreCompoundFragment::CalcEmissionProbability(const G4Fragment& fr)
50 {
51 theEmissionProbability = (Initialize(fr)) ?
52 IntegrateEmissionProbability(theMinKinEnergy, theMaxKinEnergy, fr) : 0.0;
53 /*
54 G4cout << "## G4PreCompoundFragment::CalcEmisProb "
55 << "Zf= " << fr.GetZ_asInt()
56 << " Af= " << fr.GetA_asInt()
57 << " Elow= " << theMinKinEnergy
58 << " Eup= " << theMaxKinEnergy
59 << " prob= " << theEmissionProbability
60 << " index=" << index << " Z=" << theZ << " A=" << theA
61 << G4endl;
62 */
63 return theEmissionProbability;
64 }
65
66 G4double
67 G4PreCompoundFragment::IntegrateEmissionProbability(G4double low, G4double up,
68 const G4Fragment& fr)
69 {
70 static const G4double den = 1.0/CLHEP::MeV;
71 G4double del = (up - low);
72 G4int nbins = del*den;
73 nbins = std::max(nbins, 4);
74 del /= static_cast<G4double>(nbins);
75 G4double e = low + 0.5*del;
76 probmax = ProbabilityDistributionFunction(e, fr);
77 //G4cout << " 0. e= " << e << " y= " << probmax << G4endl;
78
79 G4double sum = probmax;
80 for (G4int i=1; i<nbins; ++i) {
81 e += del;
82
83 G4double y = ProbabilityDistributionFunction(e, fr);
84 probmax = std::max(probmax, y);
85 sum += y;
86 if(y < sum*0.01) { break; }
87 //G4cout <<" "<<i<<". e= "<<e<<" y= "<<y<<" sum= "<<sum<< G4endl;
88 }
89 sum *= del;
90 //G4cout << "Evap prob: " << sum << " probmax= " << probmax << G4endl;
91 return sum;
92 }
93
94 G4double G4PreCompoundFragment::CrossSection(G4double ekin)
95 {
96 /*
97 G4cout << "G4PreCompoundFragment::CrossSection OPTxs=" << OPTxs << " E=" << ekin
98 << " resZ=" << theResZ << " resA=" << theResA << " index=" << index
99 << " fXSection:" << fXSection << G4endl;
100 */
101 // compute power once
102 if (OPTxs > 1 && 0 < index && theResA != lastA) {
103 lastA = theResA;
104 muu = G4KalbachCrossSection::ComputePowerParameter(lastA, index);
105 }
106 if (OPTxs == 0) {
107 recentXS = GetOpt0(ekin);
108 } else if (OPTxs == 1) {
109 G4int Z = std::min(theResZ, ZMAXNUCLEARDATA);
110 //G4double e = std::max(ekin, lowEnergyLimitMeV[Z]);
111 recentXS = fXSection->GetElementCrossSection(ekin, Z)/CLHEP::millibarn;
112
113 } else if (OPTxs == 2) {
114 recentXS = G4ChatterjeeCrossSection::ComputeCrossSection(ekin,
115 theCoulombBarrier,
116 theResA13, muu,
117 index, theZ, theResA);
118
119 } else {
120 recentXS = G4KalbachCrossSection::ComputeCrossSection(ekin, theCoulombBarrier,
121 theResA13, muu, index,
122 theZ, theA, theResA);
123 }
124 return recentXS;
125 }
126
127 G4double G4PreCompoundFragment::GetOpt0(G4double ekin) const
128 // OPT=0 : Dostrovski's cross section
129 {
130 G4double r0 = theParameters->GetR0()*theResA13;
131 // cross section is now given in mb (r0 is in mm) for the sake of consistency
132 // with the rest of the options
133 return 1.e+25*CLHEP::pi*r0*r0*theResA13*GetAlpha()*(1.0 + GetBeta()/ekin);
134 }
135
136 G4double G4PreCompoundFragment::SampleKineticEnergy(const G4Fragment& fragment)
137 {
138 G4double delta = theMaxKinEnergy - theMinKinEnergy;
139 static const G4double toler = 1.25;
140 probmax *= toler;
141 G4double prob, T(0.0);
142 CLHEP::HepRandomEngine* rndm = G4Random::getTheEngine();
143 G4int i;
144 for(i=0; i<100; ++i) {
145 T = theMinKinEnergy + delta*rndm->flat();
146 prob = ProbabilityDistributionFunction(T, fragment);
147 /*
148 if(prob > probmax) {
149 G4cout << "G4PreCompoundFragment WARNING: prob= " << prob
150 << " probmax= " << probmax << G4endl;
151 G4cout << "i= " << i << " Z= " << theZ << " A= " << theA
152 << " resZ= " << theResZ << " resA= " << theResA << "\n"
153 << " T= " << T << " Tmax= " << theMaxKinEnergy
154 << " Tmin= " << limit
155 << G4endl;
156 for(G4int i=0; i<N; ++i) { G4cout << " " << probability[i]; }
157 G4cout << G4endl;
158 }
159 */
160 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
161 if(probmax*rndm->flat() <= prob) { break; }
162 }
163 /*
164 G4cout << "G4PreCompoundFragment: i= " << i << " Z= " << theZ
165 << " A= " << theA <<" T(MeV)= " << T << " Emin(MeV)= "
166 << theMinKinEnergy << " Emax= " << theMaxKinEnergy << G4endl;
167 */
168 return T;
169 }
170
171