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