| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 09 Dec 2009 L Pandola First implementa
32 //
33 #include "G4PenelopeSamplingData.hh"
34
35 //....oooOO0OOooo........oooOO0OOooo........oo
36 G4PenelopeSamplingData::G4PenelopeSamplingData
37 fNP(nPoints)
38 {
39 //create vectors
40 fX = new G4DataVector();
41 fPAC = new G4DataVector();
42 fA = new G4DataVector();
43 fB = new G4DataVector();
44 fITTL = new std::vector<size_t>;
45 fITTU = new std::vector<size_t>;
46 }
47
48 //....oooOO0OOooo........oooOO0OOooo........oo
49 G4PenelopeSamplingData::~G4PenelopeSamplingDat
50 {
51 if (fX) delete fX;
52 if (fPAC) delete fPAC;
53 if (fA) delete fA;
54 if (fB) delete fB;
55 if (fITTL) delete fITTL;
56 if (fITTU) delete fITTU;
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 size_t G4PenelopeSamplingData::GetNumberOfStor
62 {
63 size_t points = fX->size();
64
65 //check everything is all right
66 if (fPAC->size() != points || fA->size() !=
67 fB->size() != points || fITTL->size() !=
68 fITTU->size() != points)
69 {
70 G4ExceptionDescription ed;
71 ed << "Data vectors look to have differe
72 G4Exception("G4PenelopeSamplingData::Get
73 FatalException,ed);
74 }
75 return points;
76 }
77
78 //....oooOO0OOooo........oooOO0OOooo........oo
79 void G4PenelopeSamplingData::Clear()
80 {
81 if (fX) delete fX;
82 if (fPAC) delete fPAC;
83 if (fA) delete fA;
84 if (fB) delete fB;
85 if (fITTL) delete fITTL;
86 if (fITTU) delete fITTU;
87 //create vectors
88 fX = new G4DataVector();
89 fPAC = new G4DataVector();
90 fA = new G4DataVector();
91 fB = new G4DataVector();
92 fITTL = new std::vector<size_t>;
93 fITTU = new std::vector<size_t>;
94 }
95
96 //....oooOO0OOooo........oooOO0OOooo........oo
97 void G4PenelopeSamplingData::AddPoint(G4double
98 size_t ITTL0,size_t ITTU0)
99 {
100 fX->push_back(x0);
101 fPAC->push_back(pac0);
102 fA->push_back(a0);
103 fB->push_back(b0);
104 fITTL->push_back(ITTL0);
105 fITTU->push_back(ITTU0);
106
107 //check how many points we do have now
108 size_t nOfPoints = GetNumberOfStoredPoints()
109
110 if (nOfPoints > ((size_t)fNP))
111 {
112 G4cout << "G4PenelopeSamplingData::AddPo
113 G4cout << "WARNING: Up to now there are
114 G4cout << "while the anticipated (declar
115 }
116 return;
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oo
120 void G4PenelopeSamplingData::DumpTable()
121 {
122
123 G4cout << "*********************************
124 G4cout << GetNumberOfStoredPoints() << " poi
125 G4cout << "*********************************
126 for (size_t i=0;i<GetNumberOfStoredPoints();
127 {
128 G4cout << i << " " << (*fX)[i] << " " <<
129 (*fB)[i] << " " << (*fITTL)[i] << " " << (*f
130 }
131 G4cout << "*********************************
132 }
133
134 //....oooOO0OOooo........oooOO0OOooo........oo
135 G4double G4PenelopeSamplingData::GetX(size_t i
136 {
137 if (index < fX->size())
138 return (*fX)[index];
139 else
140 return 0;
141 }
142
143 //....oooOO0OOooo........oooOO0OOooo........oo
144 G4double G4PenelopeSamplingData::GetPAC(size_t
145 {
146 if (index < fPAC->size())
147 return (*fPAC)[index];
148 else
149 return 0;
150 }
151
152 //....oooOO0OOooo........oooOO0OOooo........oo
153 G4double G4PenelopeSamplingData::GetA(size_t i
154 {
155 if (index < fA->size())
156 return (*fA)[index];
157 else
158 return 0;
159 }
160
161 //....oooOO0OOooo........oooOO0OOooo........oo
162 G4double G4PenelopeSamplingData::GetB(size_t i
163 {
164 if (index < fB->size())
165 return (*fB)[index];
166 else
167 return 0;
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oo
171 G4double G4PenelopeSamplingData::SampleValue(G
172 {
173 //One passes here a random number in (0,1).
174 //Notice: it possible that is between (0,b)
175 size_t points = GetNumberOfStoredPoints();
176
177 size_t itn = (size_t) (maxRand*(points-1));
178 size_t i = (*fITTL)[itn];
179 size_t j = (*fITTU)[itn];
180
181 while ((j-i) > 1)
182 {
183 size_t k = (i+j)/2;
184 if (maxRand > (*fPAC)[k])
185 i = k;
186 else
187 j = k;
188 }
189
190 //Sampling from the rational inverse cumulat
191 G4double result = 0;
192
193 G4double rr = maxRand - (*fPAC)[i];
194 if (rr > 1e-16)
195 {
196 G4double d = (*fPAC)[i+1]-(*fPAC)[i];
197 result = (*fX)[i]+
198 ((1.0+(*fA)[i]+(*fB)[i])*d*rr/
199 (d*d+((*fA)[i]*d+(*fB)[i]*rr)*rr))*((*fX)[i
200 }
201 else
202 result = (*fX)[i];
203
204 return result;
205 }
206