| 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 // Geant4 class G4HadXSHelper
28 //
29 // Author V.Ivanchenko 18.05.2022
30 //
31
32 #include "G4HadXSHelper.hh"
33 #include "G4Material.hh"
34 #include "G4MaterialTable.hh"
35 #include "G4Log.hh"
36 #include "G4Exp.hh"
37
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
39
40 static const G4double scale = 10./G4Log(10.);
41
42 std::vector<G4double>*
43 G4HadXSHelper::FindCrossSectionMax(G4HadronicProcess* p,
44 const G4ParticleDefinition* part,
45 const G4double emin,
46 const G4double emax)
47 {
48 std::vector<G4double>* ptr = nullptr;
49 if(nullptr == p || nullptr == part) { return ptr; }
50 /*
51 G4cout << "G4HadXSHelper::FindCrossSectionMax for "
52 << p->GetProcessName() << " and " << part->GetParticleName() << G4endl;
53 */
54
55 const G4MaterialTable* theMatTable = G4Material::GetMaterialTable();
56 size_t n = G4Material::GetNumberOfMaterials();
57 ptr = new std::vector<G4double>;
58 ptr->resize(n, DBL_MAX);
59
60 G4bool isPeak = false;
61 G4double ee = G4Log(emax/emin);
62 G4int nbin = G4lrint(ee*scale);
63 if(nbin < 4) { nbin = 4; }
64 G4double x = G4Exp(ee/nbin);
65
66 // first loop on existing vectors
67 for (size_t i=0; i<n; ++i) {
68 auto mat = (*theMatTable)[i];
69 G4double sm = 0.0;
70 G4double em = 0.0;
71 G4double e = emin;
72 for(G4int j=0; j<=nbin; ++j) {
73 G4double sig = p->ComputeCrossSection(part, mat, e);
74 if(sig >= sm) {
75 em = e;
76 sm = sig;
77 e = (j+1 < nbin) ? e*x : emax;
78 } else {
79 isPeak = true;
80 (*ptr)[i] = em;
81 break;
82 }
83 }
84 //G4cout << i << ". em=" << em << " sm=" << sm << G4endl;
85 }
86 // there is no peak for any couple
87 if(!isPeak) {
88 delete ptr;
89 ptr = nullptr;
90 }
91 return ptr;
92 }
93
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
95
96 std::vector<G4TwoPeaksHadXS*>*
97 G4HadXSHelper::FillPeaksStructure(G4HadronicProcess* p,
98 const G4ParticleDefinition* part,
99 const G4double emin,
100 const G4double emax)
101 {
102 std::vector<G4TwoPeaksHadXS*>* ptr = nullptr;
103 if(nullptr == p) { return ptr; }
104
105 /*
106 G4cout << "G4HadXSHelper::FillPeaksStructure for "
107 << p->GetProcessName() << " and " << part->GetParticleName() << G4endl;
108 */
109 const G4MaterialTable* theMatTable = G4Material::GetMaterialTable();
110 size_t n = G4Material::GetNumberOfMaterials();
111 ptr = new std::vector<G4TwoPeaksHadXS*>;
112 ptr->resize(n, nullptr);
113
114 G4double e1peak, e1deep, e2peak, e2deep, e3peak;
115 G4bool isDeep = false;
116
117 G4double ee = G4Log(emax/emin);
118 G4int nbin = G4lrint(ee*scale);
119 if(nbin < 4) { nbin = 4; }
120 G4double fact = G4Exp(ee/nbin);
121
122 for (size_t i=0; i<n; ++i) {
123 auto mat = (*theMatTable)[i];
124 G4double e = emin/fact;
125
126 G4double xs = 0.0;
127 e1peak = e1deep = e2peak = e2deep = e3peak = DBL_MAX;
128 for(G4int j=0; j<=nbin; ++j) {
129 e = (j+1 < nbin) ? e*fact : emax;
130 G4double ss = p->ComputeCrossSection(part, mat, e);
131 // find out 1st peak
132 if(e1peak == DBL_MAX) {
133 if(ss >= xs) {
134 xs = ss;
135 ee = e;
136 continue;
137 } else {
138 e1peak = ee;
139 }
140 }
141 // find out the deep
142 if(e1deep == DBL_MAX) {
143 if(ss <= xs) {
144 xs = ss;
145 ee = e;
146 continue;
147 } else {
148 e1deep = ee;
149 isDeep = true;
150 }
151 }
152 // find out 2nd peak
153 if(e2peak == DBL_MAX) {
154 if(ss >= xs) {
155 xs = ss;
156 ee = e;
157 continue;
158 } else {
159 e2peak = ee;
160 }
161 }
162 if(e2deep == DBL_MAX) {
163 if(ss <= xs) {
164 xs = ss;
165 ee = e;
166 continue;
167 } else {
168 e2deep = ee;
169 break;
170 }
171 }
172 // find out 3d peak
173 if(e3peak == DBL_MAX) {
174 if(ss >= xs) {
175 xs = ss;
176 ee = e;
177 continue;
178 } else {
179 e3peak = ee;
180 }
181 }
182 }
183 G4TwoPeaksHadXS* x = (*ptr)[i];
184 if(nullptr == x) {
185 x = new G4TwoPeaksHadXS();
186 (*ptr)[i] = x;
187 }
188 x->e1peak = e1peak;
189 x->e1deep = e1deep;
190 x->e2peak = e2peak;
191 x->e2deep = e2deep;
192 x->e3peak = e3peak;
193 //G4cout << "coupleIdx=" << i << " " << e1peak << " " << e1deep
194 // << " " << e2peak << " " << e2deep << " " << e3peak << G4endl;
195 }
196 // case of no 1st peak in all vectors
197 if(!isDeep) {
198 for (auto& x : *ptr) {
199 delete x;
200 }
201 delete ptr;
202 ptr = nullptr;
203 }
204 return ptr;
205 }
206
207 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
208