| 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 #include "G4AdjointCSMatrix.hh" 27 #include "G4AdjointCSMatrix.hh"
28 28
29 #include "G4AdjointInterpolator.hh" 29 #include "G4AdjointInterpolator.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 31
32 #include <iomanip> 32 #include <iomanip>
33 #include <fstream> 33 #include <fstream>
34 34
35 ////////////////////////////////////////////// 35 ///////////////////////////////////////////////////////
36 G4AdjointCSMatrix::G4AdjointCSMatrix(G4bool aB 36 G4AdjointCSMatrix::G4AdjointCSMatrix(G4bool aBool) { fScatProjToProj = aBool; }
37 37
38 ////////////////////////////////////////////// 38 ///////////////////////////////////////////////////////
39 G4AdjointCSMatrix::~G4AdjointCSMatrix() 39 G4AdjointCSMatrix::~G4AdjointCSMatrix()
40 { 40 {
41 fLogPrimEnergyVector.clear(); 41 fLogPrimEnergyVector.clear();
42 fLogCrossSectionVector.clear(); 42 fLogCrossSectionVector.clear();
43 43
44 for (auto p : fLogSecondEnergyMatrix) { 44 for (auto p : fLogSecondEnergyMatrix) {
45 p->clear(); 45 p->clear();
46 delete p; 46 delete p;
47 p = nullptr; 47 p = nullptr;
48 } 48 }
49 fLogSecondEnergyMatrix.clear(); 49 fLogSecondEnergyMatrix.clear();
50 50
51 for (auto p : fLogProbMatrix) { 51 for (auto p : fLogProbMatrix) {
52 p->clear(); 52 p->clear();
53 delete p; 53 delete p;
54 p = nullptr; 54 p = nullptr;
55 } 55 }
56 fLogProbMatrix.clear(); 56 fLogProbMatrix.clear();
57 57
58 for (auto p : fLogProbMatrixIndex) { 58 for (auto p : fLogProbMatrixIndex) {
59 if (p) { 59 if (p) {
60 p->clear(); 60 p->clear();
61 delete p; 61 delete p;
62 p = nullptr; 62 p = nullptr;
63 } 63 }
64 } 64 }
65 fLogProbMatrixIndex.clear(); 65 fLogProbMatrixIndex.clear();
66 } 66 }
67 67
68 ////////////////////////////////////////////// 68 ///////////////////////////////////////////////////////
69 void G4AdjointCSMatrix::Clear() 69 void G4AdjointCSMatrix::Clear()
70 { 70 {
71 fLogPrimEnergyVector.clear(); 71 fLogPrimEnergyVector.clear();
72 fLogCrossSectionVector.clear(); 72 fLogCrossSectionVector.clear();
73 fLogSecondEnergyMatrix.clear(); 73 fLogSecondEnergyMatrix.clear();
74 fLogProbMatrix.clear(); 74 fLogProbMatrix.clear();
75 fLogProbMatrixIndex.clear(); 75 fLogProbMatrixIndex.clear();
76 fLog0Vector.clear(); 76 fLog0Vector.clear();
77 fNbPrimEnergy = 0; 77 fNbPrimEnergy = 0;
78 } 78 }
79 79
80 ////////////////////////////////////////////// 80 ///////////////////////////////////////////////////////
81 void G4AdjointCSMatrix::AddData(G4double aLogP 81 void G4AdjointCSMatrix::AddData(G4double aLogPrimEnergy, G4double aLogCS,
82 std::vector<G4 << 82 std::vector<double>* aLogSecondEnergyVector,
83 std::vector<G4 << 83 std::vector<double>* aLogProbVector,
84 std::size_t n_ << 84 size_t n_pro_decade)
85 { 85 {
86 G4AdjointInterpolator* theInterpolator = G4A 86 G4AdjointInterpolator* theInterpolator = G4AdjointInterpolator::GetInstance();
87 87
88 // At this time we consider that the energy 88 // At this time we consider that the energy is increasing monotically
89 fLogPrimEnergyVector.push_back(aLogPrimEnerg 89 fLogPrimEnergyVector.push_back(aLogPrimEnergy);
90 fLogCrossSectionVector.push_back(aLogCS); 90 fLogCrossSectionVector.push_back(aLogCS);
91 fLogSecondEnergyMatrix.push_back(aLogSecondE 91 fLogSecondEnergyMatrix.push_back(aLogSecondEnergyVector);
92 fLogProbMatrix.push_back(aLogProbVector); 92 fLogProbMatrix.push_back(aLogProbVector);
93 93
94 std::vector<std::size_t>* aLogProbVectorInde << 94 std::vector<size_t>* aLogProbVectorIndex = nullptr;
95 95
96 if(n_pro_decade > 0 && !aLogProbVector->empt 96 if(n_pro_decade > 0 && !aLogProbVector->empty())
97 { 97 {
98 aLogProbVectorIndex = new std::vector<std: << 98 aLogProbVectorIndex = new std::vector<size_t>();
99 G4double dlog = std::log(10.) / n_pr 99 G4double dlog = std::log(10.) / n_pro_decade;
100 G4double log_val = 100 G4double log_val =
101 G4int(std::min((*aLogProbVector)[0], aLo << 101 int(std::min((*aLogProbVector)[0], aLogProbVector->back()) / dlog) * dlog;
102 fLog0Vector.push_back(log_val); 102 fLog0Vector.push_back(log_val);
103 103
104 // Loop checking, 07-Aug-2015, Vladimir Iv 104 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
105 while(log_val < 0.) 105 while(log_val < 0.)
106 { 106 {
107 aLogProbVectorIndex->push_back( 107 aLogProbVectorIndex->push_back(
108 theInterpolator->FindPosition(log_val, 108 theInterpolator->FindPosition(log_val, (*aLogProbVector)));
109 log_val += dlog; 109 log_val += dlog;
110 } 110 }
111 } 111 }
112 else 112 else
113 { 113 {
114 fLog0Vector.push_back(0.); 114 fLog0Vector.push_back(0.);
115 } 115 }
116 fLogProbMatrixIndex.push_back(aLogProbVector 116 fLogProbMatrixIndex.push_back(aLogProbVectorIndex);
117 117
118 ++fNbPrimEnergy; 118 ++fNbPrimEnergy;
119 } 119 }
120 120
121 ////////////////////////////////////////////// 121 ///////////////////////////////////////////////////////
122 G4bool G4AdjointCSMatrix::GetData(unsigned int 122 G4bool G4AdjointCSMatrix::GetData(unsigned int i, G4double& aLogPrimEnergy,
123 G4double& aL 123 G4double& aLogCS, G4double& log0,
124 std::vector< << 124 std::vector<double>*& aLogSecondEnergyVector,
125 std::vector< << 125 std::vector<double>*& aLogProbVector,
126 std::vector< << 126 std::vector<size_t>*& aLogProbVectorIndex)
127 { 127 {
128 if(i >= fNbPrimEnergy) 128 if(i >= fNbPrimEnergy)
129 return false; 129 return false;
130 aLogPrimEnergy = fLogPrimEnergyVecto 130 aLogPrimEnergy = fLogPrimEnergyVector[i];
131 aLogCS = fLogCrossSectionVec 131 aLogCS = fLogCrossSectionVector[i];
132 aLogSecondEnergyVector = fLogSecondEnergyMat 132 aLogSecondEnergyVector = fLogSecondEnergyMatrix[i];
133 aLogProbVector = fLogProbMatrix[i]; 133 aLogProbVector = fLogProbMatrix[i];
134 aLogProbVectorIndex = fLogProbMatrixIndex 134 aLogProbVectorIndex = fLogProbMatrixIndex[i];
135 log0 = fLog0Vector[i]; 135 log0 = fLog0Vector[i];
136 return true; 136 return true;
137 } 137 }
138 138
139 ////////////////////////////////////////////// 139 ///////////////////////////////////////////////////////
140 void G4AdjointCSMatrix::Write(const G4String& << 140 void G4AdjointCSMatrix::Write(G4String file_name)
141 { 141 {
142 std::fstream FileOutput(file_name, std::ios: 142 std::fstream FileOutput(file_name, std::ios::out);
143 FileOutput << std::setiosflags(std::ios::sci 143 FileOutput << std::setiosflags(std::ios::scientific);
144 FileOutput << std::setprecision(6); 144 FileOutput << std::setprecision(6);
145 FileOutput << fLogPrimEnergyVector.size() << 145 FileOutput << fLogPrimEnergyVector.size() << G4endl;
146 for(std::size_t i = 0; i < fLogPrimEnergyVec << 146 for(size_t i = 0; i < fLogPrimEnergyVector.size(); ++i)
147 { 147 {
148 FileOutput << std::exp(fLogPrimEnergyVecto 148 FileOutput << std::exp(fLogPrimEnergyVector[i]) / MeV << '\t'
149 << std::exp(fLogCrossSectionVec 149 << std::exp(fLogCrossSectionVector[i]) << G4endl;
150 std::size_t j1 = 0; << 150 size_t j1 = 0;
151 FileOutput << fLogSecondEnergyMatrix[i]->s 151 FileOutput << fLogSecondEnergyMatrix[i]->size() << G4endl;
152 for(std::size_t j = 0; j < fLogSecondEnerg << 152 for(size_t j = 0; j < fLogSecondEnergyMatrix[i]->size(); ++j)
153 { 153 {
154 FileOutput << std::exp((*fLogSecondEnerg 154 FileOutput << std::exp((*fLogSecondEnergyMatrix[i])[j]);
155 ++j1; 155 ++j1;
156 if(j1 < 10) 156 if(j1 < 10)
157 FileOutput << '\t'; 157 FileOutput << '\t';
158 else 158 else
159 { 159 {
160 FileOutput << G4endl; 160 FileOutput << G4endl;
161 j1 = 0; 161 j1 = 0;
162 } 162 }
163 } 163 }
164 if(j1 > 0) 164 if(j1 > 0)
165 FileOutput << G4endl; 165 FileOutput << G4endl;
166 j1 = 0; 166 j1 = 0;
167 FileOutput << fLogProbMatrix[i]->size() << 167 FileOutput << fLogProbMatrix[i]->size() << G4endl;
168 for(std::size_t j = 0; j < fLogProbMatrix[ << 168 for(size_t j = 0; j < fLogProbMatrix[i]->size(); ++j)
169 { 169 {
170 FileOutput << std::exp((*fLogProbMatrix[ 170 FileOutput << std::exp((*fLogProbMatrix[i])[j]);
171 ++j1; 171 ++j1;
172 if(j1 < 10) 172 if(j1 < 10)
173 FileOutput << '\t'; 173 FileOutput << '\t';
174 else 174 else
175 { 175 {
176 FileOutput << G4endl; 176 FileOutput << G4endl;
177 j1 = 0; 177 j1 = 0;
178 } 178 }
179 } 179 }
180 if(j1 > 0) 180 if(j1 > 0)
181 FileOutput << G4endl; 181 FileOutput << G4endl;
182 } 182 }
183 } 183 }
184 184
185 ////////////////////////////////////////////// 185 ///////////////////////////////////////////////////////
186 void G4AdjointCSMatrix::Read(const G4String& f << 186 void G4AdjointCSMatrix::Read(G4String file_name)
187 { 187 {
188 std::fstream FileOutput(file_name, std::ios: 188 std::fstream FileOutput(file_name, std::ios::in);
189 std::size_t n1, n2; << 189 size_t n1, n2;
190 190
191 fLogPrimEnergyVector.clear(); 191 fLogPrimEnergyVector.clear();
192 fLogCrossSectionVector.clear(); 192 fLogCrossSectionVector.clear();
193 fLogSecondEnergyMatrix.clear(); 193 fLogSecondEnergyMatrix.clear();
194 fLogProbMatrix.clear(); 194 fLogProbMatrix.clear();
195 FileOutput >> n1; 195 FileOutput >> n1;
196 for(std::size_t i = 0; i < n1; ++i) << 196 for(size_t i = 0; i < n1; ++i)
197 { 197 {
198 G4double E, CS; 198 G4double E, CS;
199 FileOutput >> E >> CS; 199 FileOutput >> E >> CS;
200 fLogPrimEnergyVector.push_back(E); 200 fLogPrimEnergyVector.push_back(E);
201 fLogCrossSectionVector.push_back(CS); 201 fLogCrossSectionVector.push_back(CS);
202 FileOutput >> n2; 202 FileOutput >> n2;
203 fLogSecondEnergyMatrix.push_back(new std:: 203 fLogSecondEnergyMatrix.push_back(new std::vector<G4double>());
204 fLogProbMatrix.push_back(new std::vector<G 204 fLogProbMatrix.push_back(new std::vector<G4double>());
205 205
206 for(std::size_t j = 0; j < n2; ++j) << 206 for(size_t j = 0; j < n2; ++j)
207 { 207 {
208 G4double E1; 208 G4double E1;
209 FileOutput >> E1; 209 FileOutput >> E1;
210 fLogSecondEnergyMatrix[i]->push_back(E1) 210 fLogSecondEnergyMatrix[i]->push_back(E1);
211 } 211 }
212 FileOutput >> n2; 212 FileOutput >> n2;
213 for(std::size_t j = 0; j < n2; ++j) << 213 for(size_t j = 0; j < n2; ++j)
214 { 214 {
215 G4double prob; 215 G4double prob;
216 FileOutput >> prob; 216 FileOutput >> prob;
217 fLogProbMatrix[i]->push_back(prob); 217 fLogProbMatrix[i]->push_back(prob);
218 } 218 }
219 } 219 }
220 } 220 }
221 221