Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/adjoint/src/G4AdjointCSMatrix.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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