00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "G4PenelopeSamplingData.hh"
00035
00036
00037 G4PenelopeSamplingData::G4PenelopeSamplingData(G4int nPoints) :
00038 np(nPoints)
00039 {
00040
00041 x = new G4DataVector();
00042 pac = new G4DataVector();
00043 a = new G4DataVector();
00044 b = new G4DataVector();
00045 ITTL = new std::vector<size_t>;
00046 ITTU = new std::vector<size_t>;
00047 }
00048
00049
00050 G4PenelopeSamplingData::~G4PenelopeSamplingData()
00051 {
00052 if (x) delete x;
00053 if (pac) delete pac;
00054 if (a) delete a;
00055 if (b) delete b;
00056 if (ITTL) delete ITTL;
00057 if (ITTU) delete ITTU;
00058 }
00059
00060
00061 size_t G4PenelopeSamplingData::GetNumberOfStoredPoints()
00062 {
00063 size_t points = x->size();
00064
00065
00066 if (pac->size() != points || a->size() != points ||
00067 b->size() != points || ITTL->size() != points ||
00068 ITTU->size() != points)
00069 {
00070 G4ExceptionDescription ed;
00071 ed << "Data vectors look to have different dimensions !" << G4endl;
00072 G4Exception("G4PenelopeSamplingData::GetNumberOfStoredPoints()","em2040",
00073 FatalException,ed);
00074 }
00075 return points;
00076 }
00077
00078
00079 void G4PenelopeSamplingData::Clear()
00080 {
00081 if (x) delete x;
00082 if (pac) delete pac;
00083 if (a) delete a;
00084 if (b) delete b;
00085 if (ITTL) delete ITTL;
00086 if (ITTU) delete ITTU;
00087
00088 x = new G4DataVector();
00089 pac = new G4DataVector();
00090 a = new G4DataVector();
00091 b = new G4DataVector();
00092 ITTL = new std::vector<size_t>;
00093 ITTU = new std::vector<size_t>;
00094 }
00095
00096
00097 void G4PenelopeSamplingData::AddPoint(G4double x0,G4double pac0,G4double a0,G4double b0,
00098 size_t ITTL0,size_t ITTU0)
00099 {
00100 x->push_back(x0);
00101 pac->push_back(pac0);
00102 a->push_back(a0);
00103 b->push_back(b0);
00104 ITTL->push_back(ITTL0);
00105 ITTU->push_back(ITTU0);
00106
00107
00108 size_t nOfPoints = GetNumberOfStoredPoints();
00109
00110 if (nOfPoints > ((size_t)np))
00111 {
00112 G4cout << "G4PenelopeSamplingData::AddPoint() " << G4endl;
00113 G4cout << "WARNING: Up to now there are " << nOfPoints << " points in the table" << G4endl;
00114 G4cout << "while the anticipated (declared) number is " << np << G4endl;
00115 }
00116 return;
00117 }
00118
00119
00120 void G4PenelopeSamplingData::DumpTable()
00121 {
00122
00123 G4cout << "*************************************************************************" << G4endl;
00124 G4cout << GetNumberOfStoredPoints() << " points" << G4endl;
00125 G4cout << "*************************************************************************" << G4endl;
00126 for (size_t i=0;i<GetNumberOfStoredPoints();i++)
00127 {
00128 G4cout << i << " " << (*x)[i] << " " << (*pac)[i] << " " << (*a)[i] << " " <<
00129 (*b)[i] << " " << (*ITTL)[i] << " " << (*ITTU)[i] << G4endl;
00130 }
00131 G4cout << "*************************************************************************" << G4endl;
00132 }
00133
00134
00135 G4double G4PenelopeSamplingData::GetX(size_t index)
00136 {
00137 if (index < x->size())
00138 return (*x)[index];
00139 else
00140 return 0;
00141 }
00142
00143
00144 G4double G4PenelopeSamplingData::GetPAC(size_t index)
00145 {
00146 if (index < pac->size())
00147 return (*pac)[index];
00148 else
00149 return 0;
00150 }
00151
00152
00153 G4double G4PenelopeSamplingData::GetA(size_t index)
00154 {
00155 if (index < a->size())
00156 return (*a)[index];
00157 else
00158 return 0;
00159 }
00160
00161
00162 G4double G4PenelopeSamplingData::GetB(size_t index)
00163 {
00164 if (index < b->size())
00165 return (*b)[index];
00166 else
00167 return 0;
00168 }
00169
00170
00171 G4double G4PenelopeSamplingData::SampleValue(G4double maxRand)
00172 {
00173
00174
00175 size_t points = GetNumberOfStoredPoints();
00176
00177 size_t itn = (size_t) (maxRand*(points-1));
00178 size_t i = (*ITTL)[itn];
00179 size_t j = (*ITTU)[itn];
00180
00181 while ((j-i) > 1)
00182 {
00183 size_t k = (i+j)/2;
00184 if (maxRand > (*pac)[k])
00185 i = k;
00186 else
00187 j = k;
00188 }
00189
00190
00191 G4double result = 0;
00192
00193 G4double rr = maxRand - (*pac)[i];
00194 if (rr > 1e-16)
00195 {
00196 G4double d = (*pac)[i+1]-(*pac)[i];
00197 result = (*x)[i]+
00198 ((1.0+(*a)[i]+(*b)[i])*d*rr/
00199 (d*d+((*a)[i]*d+(*b)[i]*rr)*rr))*((*x)[i+1]-(*x)[i]);
00200 }
00201 else
00202 result = (*x)[i];
00203
00204 return result;
00205 }