| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Geant4 class G4HadXSHelper 27 // Geant4 class G4HadXSHelper
28 // 28 //
29 // Author V.Ivanchenko 18.05.2022 29 // Author V.Ivanchenko 18.05.2022
30 // 30 //
31 31
32 #include "G4HadXSHelper.hh" 32 #include "G4HadXSHelper.hh"
33 #include "G4Material.hh" 33 #include "G4Material.hh"
34 #include "G4MaterialTable.hh" 34 #include "G4MaterialTable.hh"
35 #include "G4Log.hh" 35 #include "G4Log.hh"
36 #include "G4Exp.hh" 36 #include "G4Exp.hh"
37 37
38 //....oooOO0OOooo........oooOO0OOooo........oo 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
39 39
40 static const G4double scale = 10./G4Log(10.); 40 static const G4double scale = 10./G4Log(10.);
41 41
42 std::vector<G4double>* 42 std::vector<G4double>*
43 G4HadXSHelper::FindCrossSectionMax(G4HadronicP 43 G4HadXSHelper::FindCrossSectionMax(G4HadronicProcess* p,
44 const G4Par 44 const G4ParticleDefinition* part,
45 const G4dou 45 const G4double emin,
46 const G4dou 46 const G4double emax)
47 { 47 {
48 std::vector<G4double>* ptr = nullptr; 48 std::vector<G4double>* ptr = nullptr;
49 if(nullptr == p || nullptr == part) { return 49 if(nullptr == p || nullptr == part) { return ptr; }
50 /* 50 /*
51 G4cout << "G4HadXSHelper::FindCrossSectionMa 51 G4cout << "G4HadXSHelper::FindCrossSectionMax for "
52 << p->GetProcessName() << " and " << part-> 52 << p->GetProcessName() << " and " << part->GetParticleName() << G4endl;
53 */ 53 */
54 54
55 const G4MaterialTable* theMatTable = G4Mater 55 const G4MaterialTable* theMatTable = G4Material::GetMaterialTable();
56 size_t n = G4Material::GetNumberOfMaterials( 56 size_t n = G4Material::GetNumberOfMaterials();
57 ptr = new std::vector<G4double>; 57 ptr = new std::vector<G4double>;
58 ptr->resize(n, DBL_MAX); 58 ptr->resize(n, DBL_MAX);
59 59
60 G4bool isPeak = false; 60 G4bool isPeak = false;
61 G4double ee = G4Log(emax/emin); 61 G4double ee = G4Log(emax/emin);
62 G4int nbin = G4lrint(ee*scale); 62 G4int nbin = G4lrint(ee*scale);
63 if(nbin < 4) { nbin = 4; } 63 if(nbin < 4) { nbin = 4; }
64 G4double x = G4Exp(ee/nbin); 64 G4double x = G4Exp(ee/nbin);
65 65
66 // first loop on existing vectors 66 // first loop on existing vectors
67 for (size_t i=0; i<n; ++i) { 67 for (size_t i=0; i<n; ++i) {
68 auto mat = (*theMatTable)[i]; 68 auto mat = (*theMatTable)[i];
69 G4double sm = 0.0; 69 G4double sm = 0.0;
70 G4double em = 0.0; 70 G4double em = 0.0;
71 G4double e = emin; 71 G4double e = emin;
72 for(G4int j=0; j<=nbin; ++j) { 72 for(G4int j=0; j<=nbin; ++j) {
73 G4double sig = p->ComputeCrossSection(pa 73 G4double sig = p->ComputeCrossSection(part, mat, e);
74 if(sig >= sm) { 74 if(sig >= sm) {
75 em = e; 75 em = e;
76 sm = sig; 76 sm = sig;
77 e = (j+1 < nbin) ? e*x : emax; 77 e = (j+1 < nbin) ? e*x : emax;
78 } else { 78 } else {
79 isPeak = true; 79 isPeak = true;
80 (*ptr)[i] = em; 80 (*ptr)[i] = em;
81 break; 81 break;
82 } 82 }
83 } 83 }
84 //G4cout << i << ". em=" << em << " sm=" 84 //G4cout << i << ". em=" << em << " sm=" << sm << G4endl;
85 } 85 }
86 // there is no peak for any couple 86 // there is no peak for any couple
87 if(!isPeak) { 87 if(!isPeak) {
88 delete ptr; 88 delete ptr;
89 ptr = nullptr; 89 ptr = nullptr;
90 } 90 }
91 return ptr; 91 return ptr;
92 } 92 }
93 93
94 //....oooOO0OOooo........oooOO0OOooo........oo 94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
95 95
96 std::vector<G4TwoPeaksHadXS*>* 96 std::vector<G4TwoPeaksHadXS*>*
97 G4HadXSHelper::FillPeaksStructure(G4HadronicPr 97 G4HadXSHelper::FillPeaksStructure(G4HadronicProcess* p,
98 const G4ParticleDefinition* part, 98 const G4ParticleDefinition* part,
99 const G4double emin, 99 const G4double emin,
100 const G4double emax) 100 const G4double emax)
101 { 101 {
102 std::vector<G4TwoPeaksHadXS*>* ptr = nullptr 102 std::vector<G4TwoPeaksHadXS*>* ptr = nullptr;
103 if(nullptr == p) { return ptr; } 103 if(nullptr == p) { return ptr; }
104 104
105 /* 105 /*
106 G4cout << "G4HadXSHelper::FillPeaksStructure 106 G4cout << "G4HadXSHelper::FillPeaksStructure for "
107 << p->GetProcessName() << " and " << part-> 107 << p->GetProcessName() << " and " << part->GetParticleName() << G4endl;
108 */ 108 */
109 const G4MaterialTable* theMatTable = G4Mater 109 const G4MaterialTable* theMatTable = G4Material::GetMaterialTable();
110 size_t n = G4Material::GetNumberOfMaterials( 110 size_t n = G4Material::GetNumberOfMaterials();
111 ptr = new std::vector<G4TwoPeaksHadXS*>; 111 ptr = new std::vector<G4TwoPeaksHadXS*>;
112 ptr->resize(n, nullptr); 112 ptr->resize(n, nullptr);
113 113
114 G4double e1peak, e1deep, e2peak, e2deep, e3p 114 G4double e1peak, e1deep, e2peak, e2deep, e3peak;
115 G4bool isDeep = false; 115 G4bool isDeep = false;
116 116
117 G4double ee = G4Log(emax/emin); 117 G4double ee = G4Log(emax/emin);
118 G4int nbin = G4lrint(ee*scale); 118 G4int nbin = G4lrint(ee*scale);
119 if(nbin < 4) { nbin = 4; } 119 if(nbin < 4) { nbin = 4; }
120 G4double fact = G4Exp(ee/nbin); 120 G4double fact = G4Exp(ee/nbin);
121 121
122 for (size_t i=0; i<n; ++i) { 122 for (size_t i=0; i<n; ++i) {
123 auto mat = (*theMatTable)[i]; 123 auto mat = (*theMatTable)[i];
124 G4double e = emin/fact; 124 G4double e = emin/fact;
125 125
126 G4double xs = 0.0; 126 G4double xs = 0.0;
127 e1peak = e1deep = e2peak = e2deep = e3peak 127 e1peak = e1deep = e2peak = e2deep = e3peak = DBL_MAX;
128 for(G4int j=0; j<=nbin; ++j) { 128 for(G4int j=0; j<=nbin; ++j) {
129 e = (j+1 < nbin) ? e*fact : emax; 129 e = (j+1 < nbin) ? e*fact : emax;
130 G4double ss = p->ComputeCrossSection(par 130 G4double ss = p->ComputeCrossSection(part, mat, e);
131 // find out 1st peak 131 // find out 1st peak
132 if(e1peak == DBL_MAX) { 132 if(e1peak == DBL_MAX) {
133 if(ss >= xs) { 133 if(ss >= xs) {
134 xs = ss; 134 xs = ss;
135 ee = e; 135 ee = e;
136 continue; 136 continue;
137 } else { 137 } else {
138 e1peak = ee; 138 e1peak = ee;
139 } 139 }
140 } 140 }
141 // find out the deep 141 // find out the deep
142 if(e1deep == DBL_MAX) { 142 if(e1deep == DBL_MAX) {
143 if(ss <= xs) { 143 if(ss <= xs) {
144 xs = ss; 144 xs = ss;
145 ee = e; 145 ee = e;
146 continue; 146 continue;
147 } else { 147 } else {
148 e1deep = ee; 148 e1deep = ee;
149 isDeep = true; 149 isDeep = true;
150 } 150 }
151 } 151 }
152 // find out 2nd peak 152 // find out 2nd peak
153 if(e2peak == DBL_MAX) { 153 if(e2peak == DBL_MAX) {
154 if(ss >= xs) { 154 if(ss >= xs) {
155 xs = ss; 155 xs = ss;
156 ee = e; 156 ee = e;
157 continue; 157 continue;
158 } else { 158 } else {
159 e2peak = ee; 159 e2peak = ee;
160 } 160 }
161 } 161 }
162 if(e2deep == DBL_MAX) { 162 if(e2deep == DBL_MAX) {
163 if(ss <= xs) { 163 if(ss <= xs) {
164 xs = ss; 164 xs = ss;
165 ee = e; 165 ee = e;
166 continue; 166 continue;
167 } else { 167 } else {
168 e2deep = ee; 168 e2deep = ee;
169 break; 169 break;
170 } 170 }
171 } 171 }
172 // find out 3d peak 172 // find out 3d peak
173 if(e3peak == DBL_MAX) { 173 if(e3peak == DBL_MAX) {
174 if(ss >= xs) { 174 if(ss >= xs) {
175 xs = ss; 175 xs = ss;
176 ee = e; 176 ee = e;
177 continue; 177 continue;
178 } else { 178 } else {
179 e3peak = ee; 179 e3peak = ee;
180 } 180 }
181 } 181 }
182 } 182 }
183 G4TwoPeaksHadXS* x = (*ptr)[i]; 183 G4TwoPeaksHadXS* x = (*ptr)[i];
184 if(nullptr == x) { 184 if(nullptr == x) {
185 x = new G4TwoPeaksHadXS(); 185 x = new G4TwoPeaksHadXS();
186 (*ptr)[i] = x; 186 (*ptr)[i] = x;
187 } 187 }
188 x->e1peak = e1peak; 188 x->e1peak = e1peak;
189 x->e1deep = e1deep; 189 x->e1deep = e1deep;
190 x->e2peak = e2peak; 190 x->e2peak = e2peak;
191 x->e2deep = e2deep; 191 x->e2deep = e2deep;
192 x->e3peak = e3peak; 192 x->e3peak = e3peak;
193 //G4cout << "coupleIdx=" << i << " " << e1 193 //G4cout << "coupleIdx=" << i << " " << e1peak << " " << e1deep
194 // << " " << e2peak << " " << e2deep << " 194 // << " " << e2peak << " " << e2deep << " " << e3peak << G4endl;
195 } 195 }
196 // case of no 1st peak in all vectors 196 // case of no 1st peak in all vectors
197 if(!isDeep) { 197 if(!isDeep) {
198 for (auto& x : *ptr) { <<
199 delete x; <<
200 } <<
201 delete ptr; 198 delete ptr;
202 ptr = nullptr; 199 ptr = nullptr;
203 } 200 }
204 return ptr; 201 return ptr;
205 } 202 }
206 203
207 //....oooOO0OOooo........oooOO0OOooo........oo 204 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
208 205