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